ПРИМЕНЕНИЕ МЕТОДОВ ЧАСТИЦ ПРИ БАЙЕСОВСКОМ ОЦЕНИВАНИИ ВРЕМЕННЫХ РЯДОВ ЧИСЛЕННОСТИ И УЛОВОВ ПО ВОЗРАСТАМ

ВОПРОСЫ РЫБОЛОВСТВА, 2020. Том 21. №2. С. 235-249 PROBLEMS OF FISHIRIES, 2020. Vol. 21. №2. P. 235-249

ДИНАМИКА ЧИСЛЕННОСТИ

УДК 519.9

ПРИМЕНЕНИЕ МЕТОДОВ ЧАСТИЦ ПРИ БАЙЕСОВСКОМ ОЦЕНИВАНИИ ВРЕМЕННЫХ РЯДОВ ЧИСЛЕННОСТИ И УЛОВОВ ПО ВОЗРАСТАМ

© 2020 г. И. И. Шевченко

Тихоокеанский филиал Всероссийского научно-исследовательского института рыбного хозяйства и океанографии (ТИНРО), Владивосток, 690091 E-mail: igor.shevchenko@tinro-center.ru

Поступила в редакцию 14.10.2019 г.

Добавление случайных аддитивных ошибок в уравнения когорт и уловов, которые записаны в логарифмической форме, позволяют представить динамику логарифмов численности (ненаблюдаемых состояний) в виде скрытых марковских процессов первого порядка, где в качестве наблюдений выступают логарифмы уловов. Для таких моделей целесообразно рассматривать задачи байесовского (оптимального) оценивания, которые заключаются в построении функций плотности распределения вероятностей состояний при полученных наблюдениях. Решение уравнений, которым удовлетворяют искомые плотности, в аналитическом виде возможно только для ограниченного класса моделей (линейных гауссовских). При использовании методов Монте-Карло замкнутая форма вычислений заменяется на формирование выборок с соответствующими распределениями и их статистическую обработку. Поскольку получать напрямую выборки анализируемых случайных величин практически невозможно, используются приближенные последовательные вычисления и вспомогательные распределения. В работе описываются некоторые алгоритмы байесовского оценивания временных рядов численности и уловов по возрастам с последовательным применением методов частиц. Пакет Fishmetica пополнен соответствующими функциями. Влияние количества частиц на параметры результирующих апостериорных плотностей проиллюстрированы расчетами в среде Julia для тестового набора данных при различном количестве используемых частиц.

ВВЕДЕНИЕ

При оценке промысловых запасов модели в пространстве состояний используются для того, чтобы учесть ошибки обработки и наблюдений при анализе доступных временных рядов. В типичных информационных ситуациях модели строятся, например, по рядам уловов и учетных съемок (индикаторов численности). Параметры модели и состояния могут оцениваться одновременно в результате максимизации критерия правдоподобия, а при построении моделей, как правило, предполагается нормальность невязок (Schnute, 1994; Valpine,

Hilborn, 2005; Aldrin et al., 2019; Nielsen, Berg, 2019).

Автором (Шевченко, 2017; 2019) описана методика для сопоставления когортам скрытых марковских моделей по накопленным временным рядам численности и уловов по возрастам. Для этого в уравнения когорты и уравнения Баранова, представленные в логарифмической форме, добавляются аддитивные ошибки, которые имеют плотности распределения Гаусса или Лапласа. Параметры зависимостей, которые описывают коэффициенты естественной и промысловой смертности, и параметры плотностей распределения невязок рассчитываются с использованием метода максимального правдоподобия. В качестве состояний выступают логарифмы численности ^a, а в качестве наблюдений — логарифмы уловов rf в возрасте a , a е ab : ae,

где a — минимальный, а a — максимальный возраст особей, входящих в когорту. Модифицированные уравнения когорты задают стохастические динамические процессы, а уравнения Баранова — правдоподобие наблюдений. В рамках байесовского подхода вся информация о скрытых состояниях при полученных наблюдениях rf определяется апостериорными плотностями вида

.....П ), ax, ay е ab : ae. В задачах

фильтрации (реконструирования состояний по наблюдениям, доступным на текущий момент), предсказания (предварительная оценка состояний по наблюдениям, доступным на текущий момент) и сглаживания (реконструирования состояний по всем доступным наблюдениям) требуется оценивать эти плотности при определенных значениями ax и ay .

За исключением специальных случаев (линейных гауссовских моделей) апостериорные байесовские распределения не могут быть оценены аналитически. Поэтому с средины 1960-х годов предпринимались попытки построения аппроксимаций для таких распределений. В 1970-х годах разрабатывались разнообразные алгоритмы с применением методов Монте-Карло, но они не были востребованы до момента появления высокопроизводительной вычислительной техники в 1980-х годах. В последующий период основные усилия направлялись на развитие приближенных методов для случаев нелинейных и/или негауссов-ских моделей (Gordon et al., 1993; Kitagawa, 1996; Hürzeler, Künsch, 1998; Russell, Norvig, 2010; Godsill et al., 2004; Särkkä, 2013).

В работе приводится алгоритмы, предназначенные для проведения последовательного численного оценивания апостериорных плотностей при фильтрации, предсказании и сглаживании с использованием методов частиц. Пакет Fishmetica был дополнен соответствующими функциями. Приведены примеры расчетов в среде Julia байесовских оценок апостериорных распределений с использованием временных рядов из первого тестового набора данных из ВНИРО для программ, применяемых в институтах отрасли при моделировании динамики запасов. Построены графики байесовских оценок логарифмов численности, значений параметров, гистограмм и плотностей распределений. Продемонстрировано влияние количества используемых частиц на параметры результирующих апостериорных плотностей.

МАТЕРИАЛ И МЕТОДИКА

Рассматриваются промысловые запасы, для которых проведено оценивание численности и уловов за некоторое количество последовательных промысловых сезонов. Для извлечения информации из соответствующих временных рядов возможно использование моделей, которые базируются на стохастических вариантах уравнений когорт и уловов. Для таких моделей описана методика определения параметров с применением критериев максимального правдоподобия при некоторых предположениях о виде зависимостей коэффициентов естественной и промысловой смертности, а также о типах распределений невязок (Шевченко, 2017). В статье И. И. Шевченко (2019) для моделей когорт в виде байесовских сетей приведены уравнения, которым удовлетворяют апостериорные плотности при фильтрации, предсказании и сглаживании. Для случая линейных гауссовских моделей переходов и наблюдений для этих уравнений можно получить решения в аналитической форме. Для того чтобы проводить байесовское оценивание в общем случае, последовательно применяется метод Монте-Карло: формируются выборки, которые имеют распределения, приближающие целевые, и проводится их статистическая обработка. В данном разделе приводятся основные предположения, описываются методика параметрической идентификации и некоторые алгоритмы последовательного байесовского оценивания временных рядов численности и уловов по возрастам с применением методов частиц.

Параметрическая идентификация

Пусть xj и yj — заданные оценки численности и выловов в год t в зависимости от возраста рыб j, где численность измеряется в штуках, а уловы в весовых единицах, и wj — коэффициент для пересчета весовых единиц измерений в штуки, tmax — число последовательных лет, для которых доступны данные наблюдений, jmx — принятый предельный учетный возраст популяции. Через mj и fj обозначаются коэффициенты естественной и промысловой смертности, а через zj и gj — производные коэффициенты полной смертности и пересчета уловов, который определяет долю возрастной группы в улове. При использовании уравнений когорты и уловов в качестве описания механизмов взаимосвязи численности и уловов модель полностью определяется матрицами коэффициентов m и f . Для их оценивания можно использовать имеющиеся временные ряды, которые представляются с использованием матриц x, y и w. Однако количество элементов неизвестных и доступных матриц практически совпадает. Для того чтобы избавиться от переопределенности и сократить количество неизвестных параметров, считается, что каждый коэффициент естественной смертности mj представим в виде произведения когортного фактора hj, который является оценкой неблагоприятности условий выживания когорты возраста j в год t, и составляющей, зависящей исключительно от возраста mj, который задается в виде U-об-разной функции. Коэффициент промысловой смертности ftj вычисляется как произведение возрастной селективности лова sj, которая описывается некоторой сигмоидой, и доли f, определяющей скорость убыли численности (интенсивность промысла). Ряд наблюдений динамики промыслового запаса за t

" Г max

лет с максимальным возрастом jmax содержит данные о (tmax + jmax -1) когортах. С каждой такой когортой связывается физический год yr, который соответствует возрасту в один год входящих в нее рыб. Через ayr и aeyr обозначаются начальные и конечные возрастные классы соответствующей когорты C . Она

описывается заданными значениями численности Х"г, уловов У"г, весов , и с ней связаны рассчитываемые коэффициенты естественной Ма и промысловой Р смертности,

уг Г уг Г >

а также коэффициенты полной смертности

г и пересчета уловов Ьуг для возрастного класса, а е аь : ае . Имеется взаимно одног г

значное соответствие между исходным (матричными: х&,у& и т.д.) и когортными (векторными: Хуг,Ууаг и т.д.) представлениями.

Обозначим через <Еа и уа случайные

!г уг ^ уг ^

величины, определяющие значения логарифмов численности и уловов соответственно. Если граничные условия заданы при а — аЬг _1, динамику и П, можно представить в виде байесовской сети, структура которой изображена на рис. 1. Стохастические варианты уравнения когорты

^+1 = ^ _ ^(в) + с(в) (1)

& уг & уг уг \\ т / \\ X / V &

определяют модели перехода, а уравнения Баранова

п — Г _¡П^1 _ ¡пСа (в ) + у(в ) (2)

&уг ^уг уг уг \\ т/ \\ у) У >

модели наблюдения, где вектор вт включает неизвестные параметры модели И&, т&, в& и /1, а вх и ву — это векторы параметров функций плотностей распределений случайных величин и .

В предположении о лапласовском или гауссовском распределениях с и V оптимальные векторы параметров модели динамики в_ —6_ и параметров распределений

в — 0 и ву — ву могут быть оценены с применением метода максимального правдоподобия, (Шевченко, 2017). Будем в дальнейшем использовать следующие обозначения:

¿а — га(в ), 6а — с(в ), с — с(в),

уг угу т)7 уг угу ту \\ х/

(о. )•

a е a : a .

Последовательное байесовское оценивание с использованием методов частиц

Пусть yr1 — первый физический год наблюдения, yr0 = yr1 - jmax +1, yn = yr +1 — j +1 yn = yr +1 — 1, (Шев•7 2 1 max ^max 3 1 max 7 V

ченко, 2017). Зафиксируем исследуемую когорту C r, yr е yr0 : yr3, и будем опускать индекс уг во всех обозначениях, относящихся к ней. Функция плотности случайной величины £ обозначается как Р (Г) или РГ. Значение РГ (Г = х) функции плотности вероятности некоторой случайной величины £ в точке х обозначается как Р (Г = х) или Р (х), если очевидно, что речь идет именно о Гд . В аналогичных ситуациях применяются упрощенные обозначения для условных вероятностей случайных величин (при задании значений других случайных величин), например:

Р(ха | уа) = Р(Г = ха= уа), Р(ха | у-) = Р(Г = ха П = у.....п = уа). (3)

Уравнения перехода и наблюдений (1) и (2) с подставленными значениями 0т = 0т, дх = 0х и 0 у = 0 , которые получены с применением метода максимального правдоподобия, в общем случае записываются как

Г+1 = г (Г,6), а е а : а _1, (4) п = ь.а (Г V), а е аь : ае, (5)

где Га : Я х Я ^ Я, ка : Я х Я ^ Я ,0,1/- независимые от Г,г/ и друг от друга случайные величины с известными плотностями распределения Р (6 ) и Р (1/). Плотность начального распределения обозначается как Р(г -1).

Для любой выбранной когорты соответствующие уравнения (1) и (2) после подстановок записываются как

Г+1 = д _ ^а +6/, (6)

Па = Г _¡г№а _ ¡пба +1/. (7)

Стохастические модели переходов и наблюдений (4) и (5) задают марковские последовательности с условными плотностями

Р(^а+1 | Г) = = \\8[^+х _Га (Г, wa))Р(юа)¿юа, (8)

Р(п IГ) =

= 1з{па_1аа(Гпа))Р^а)¿и , (9)

где £(•) — функция Дирака1, поскольку, например,

Р(Га+11Г) = = \\ Р(Га+11 Г,юа)Р(ша I д №, (10)

а £а £а+1

значение ю не зависит от д , а значение д при известных Г,юа вычисляется по (4).

Обратная статистическая задача состоит в оценивании скрытых состояний при условии, что заданы результаты наблюдений. В байесовском контексте это означает вычисление совместных апостериорных распределений всех состояний при полученных измерениях, например, оценивание распределений Р(Г _1а | уа а), а е аь _ 1: ае. Недостаток этого подхода заключается в том, что когда добавляются новые наблюдения, оценку треИспользуется фильтрующее свойство -функции: f (х) (х-х0) dx=f (х0).

буется полностью пересчитывать. Поэтому Пусть [га1 }е1:ЛГ — выборка из N > 0

приходится отказываться от полного апосте- значений случайной величины, которая

риорного оценивания в пользу оценки более имеет плотность распределения, аппроксипростых распределений. мирующую Р(Г | уа а). Рассмотрим задачу

Рассмотрим три задачи оценива- генерирования выборки {гО^1^ из N значения распределений численности Г при ний случайной величины Г I у" , которая известных значениях уловов г/а а = уа а, имеет плотность распределения, аппрокси-где ах = ау = аг е аь _ 1: ае для фильтрации, мирующую Р(Г+1 | уа :а+1). Ее можно решать,

ax = a + k,a < jmax,k e 1: jmax — a ,ay = ae для например, с использованием повторных выпредсказания и ax = as e ab — 1: ae ,ay = ae для борок на основе значимости, когда в качестве

сглаживания. предварительного распределения выступает

одношаговое предсказание (12), а в качеФильтрация стве инструментального — распределение

При фильтрации оценивается плот- наблюдения (9) (Sarkka, 2013). Процедура

ность распределения lga | ya a, a e ab — 1: ae. генерации начинается при a = ab — 1 с выборДля того чтобы сформировать рекуррентные ки {z° — }ie1N по распределению р —1)

(последовательные) соотношения для плот- (a + 1)-ом шаге, (a +1 )e ab : af, для каждоности распределения случайной величины го i e N по z" сначала генерируется элемент

Sa+1\\ ab a+1 a+1

g | y , можно воспользоваться выраже- z через уравнение состояния так, что если

нием с использованием плотности распределения

b P(S ) выбрано значение wa, то (см. (4))

П/ a+1 \\ £a+1\\0/Sa+1 I ab :a \\ V ) * t > \\ V //

P(r+1\\ / a+1) = P(y 1 g )P(g 1 У ) (11) a+1 ra (a a \\ „„ ,

P(g 1 У ) = P(ya+1| yabaa ) , (11) £ = f(z i-Wi)" (14)

Ясно, что величины Zi+1 с правдопогде плотности y"+1 | определяются моде- добием la+1 = P(zia1 \\ za) = P((fa) 1( zta+1 \\ za))

лями наблюдения, а плотности одношагово- независимо распределены с плотностью

го предсказания lga+x | ya a и знаменатель (11) P(ga+i | ya a) . Затем по модели наблюдения

рассчитываются как маргинальные (остаточ- (9) при известном П+1 = У +1 оценивается

ные) плотности: правдоподобие каждого полученного значения z/+1 по формуле (Gordon, Salmond, 1993)

P(g | ya :a)=

Q/ a+1 I ab :a\\

= J P(gx\\xa)P(X" | y"h)dxa, (12) 1Г1 = P(z )t (15)

Z ,JP( zja+11 ya:a)

P(ya+V:a) = j

= J P(ya+1 | xa+1)P(xa+1 | yaba)dxa+1. (13).

Выборки случайных величин могут li zN zP(y a+1 | z a+1)P(za+11 "b&"4 , (16)

ювятьгя с использованием их плот- j j

или (Doucet et al., 2000)

п/ а+1 I а+1\\Л/ а+1 I а :а\\

= _ Р(у | г. )Р(г. | у )

генерироваться с использованием их плотностей распределения, а по выборке можно

приближенно восстановить плотность рас- где Р(уа+1 | д^1) = Р((Ь.а+1) 1 (уа+1,1)) ,

пределения (как гистограмму). Поскольку I е , (Ь.а+ )_ — обратная к На функция,

аналитические вычисления плотностей апо- которая возвращает значение V по заданным

стериорных вероятностей возможны только значениям Г и П ; см. (5). На заклюдля специальных случаев, используются ре- чительном этапе генерируется выборка слукуррентные процедуры, которые по выбор- чайной дискретной величины, принимающей

кам априорных позволяет получать выборки значение с правдоподобием . При

апостериорных величин. N плотность ее распределения приближает плотность распределения случайной Выборка на (к +1) -шаге генерируетвеличины | уа а+1: ся с использованием соотношения (см. (4)):

Г)/ а+1 I аЬ :а+1\\ а+1 а+1 \\ / ЛЧ\\ - а&+к+1_ га&+к ( а&+к а& +к\\ /ПЛ\\

Р(х 1 у )«X,&=¡1 х _- ). (17) - = I (- ,Ю.1 (21)

Описанный метод имеет вычисли- где юа1+к распределена с плотностью Р(с),

тельную сложность О ((ае _ аь) N). а пересчет правдоподобия осуществляется по

Оценивание точности описанного ал- формуле: горитма для конечных значений N представляет собой самостоятельную задачу (Спэап, а,+к1 Р( га&+к+1 | г"&+к) I"&+к

Doucet, 2019). Для некоторых приложений иN р(+к+1 i z a&+к) ^ a&

на практике применялись значения N поряд- ¿-ч=1 ^

ка 1000-10000 (Kitagawa, 1996). В процессе

выполнения алгоритма может воз ни кат ь про- i е 1: N. (22)

блема вырождения, когда почти все частицы

имеют близкие к нулю вероятности. Она раз- Генерация по (21) и (22) начинаются

решается проведением повторной выборки. с к = 0, когда правдоподобие {Х +1} ie1N для Оценка числа эффективных частиц может P(£" 1 | ya" ) определяется соотношением

проводится с использованием формулы (14) для одношагового предсказания.

Другой вариант генерации может

п" =ф_(18) быть осуществлен через процедуру интегри(Г. )2 рования с использованием метода МонтеКарло. Пусть {z_a +к} 1N — выборка из N

ia « . « « «

где I. — нормализованный вес частицы i значений случайной величины, которая име-шаге a. Повторная выборка осуществляется ет плотность распределения, аппроксимив случаях, когда п" значительно меньше N, рующую P(£" +к | y" " ), с весами {I." +к} ,е1-.

например, п" < N /10 (Smith, Gelfand, 1992; Нетрудно видеть, что с использованием (19) Gordon et al., 1993; Kitagawa, 1996; Doucet веса для {z/ +k+1} ie1N можно пересчитывать

et al., 2000; Sarkka, 2013). по формуле:

&+k+1

Предсказание

Предсказание представляет собой ^ (Р-.а&+к+1\\ га&+к) Iа&+к (га&+к _ га&+к)

фильтрацию без добавления новых наблю - —-----т----г, (23)

„ „ Р а& +к+1 | а& +к ) ¡а& +к / а& +к _ а& +к \\

дений и заключается в оценке распределений ^&- 1 -п- Щ ^-+1 -п- } состояний +к | уа &а для к е 1: &тах _ а& при

&тах > а&. Процедура для предсказания бази- где {—а +к}-еШ есть упорядоченная по возраруется на рекуррентном соотношении: станию последовательность {—а +к}

Р(£а + + I ya :a ) = Сглаживание

Под сглаживанием мы понимаем = | | +k )Р(/+k | yab:a& )d/+k. (19) оценку распределения | yaV , a e ab -1: a. Процедура для сглаживания может ба-Пусть приближенное (эмпириче- зироваться на различных рекуррентных ское) представление плотности одношагово- соотношениях. В работе И. И. Шевченко го предсказывания на k -шаге, k e 1: jmax - ae, (2019) был рассмотрен случай представле-имеет вид: ния сглаживания через прямую и обратную

Р( a- +k I ab :a e ) ^y N . + ,!, + , __ + , ) фильтраЦию. АнаЛогично можно получить

(x 1 y ) ~yj=t-j \\x Zj )&( ) соотношение:

P(l+1 | y"") х| y""" ) дится к вычислению новых весов для частиц,

полученных при фильтрации рассмотренным P(l+1 | y"" ")P(ya+1:"e |l+1), (24) выше методом последовательной повторной

выборки с учетом значимости. где P(Ç"+ | y" " ) представляет собой плот- Пусть вычислены пары l" и z",

ность одношагового предсказания, а а G а" - 1 : a&, i G 1:N, определяющие прибли-p(y"+i a | ) вычисляется в обратном направ- женную плотность отфильтрованного солении (от k — ае -1 до k — а) с использовани- стояния. Определим новые веса l"а для —а, ем следующих зависимостей: а G а" -1 : a&, i G 1 : N для приближенной плотР(у"& | Е"& ) = Р(у"& \\ i ) Р(ук"& | i ) — ностисглаженногосостоянияследующимобразом. Для а — а& полагаем

— P(yk )P(yk+1:"& |lk), P(yk+1:a" |lk) — i =lft (30)

= | P(yk+1:"& | xk+1)P(xk+1 | i)dxk+1. (25) а затем от а — а& -1 до а — а" вычисляем Воспользуемся маргинальным соот- l"P(z "+1 | z")

l""& — у i"+1|a&—i_}=i—i G 1: n (31)

ием: i & У ГР(za+1|z")& & &

¿—¡m m 1 —m &

Р(Е" | y" " ) — J P(l | x"+1, y")dx"+1. (26) Тогда получаем, что

Поскольку i" не зависит от y"1" при Р(x | y )«У¡j" x - — )&

заданном i , с применением формулы Бай- а G аъ -1 : а & (32) еса получаем:

Р(£"+1 | £")Р(£" | "" ") Для того чтобы обосновать описанР(£" 11"1,y" " ) —-|--. (27) ный алгоритм, отметим, что в силу соотноО/р"+1 I " « /о\\ /&fix

P(l | y ) шений (о) и (17) выполняется приближенное

Перепишем теперь плотность со- равенство:

вместного распределения и i при задана" :а

ном y :

ношением:

p(xa+1| ya:a ) «у n—1lmp(xa+1| xm ). (33)

P(E,a | y"b:a&) = Крометого,

= a ) p( +, | ba ) (28) f )P(*°+1|y- ) и

n/£"+ll ab :a\\ ^ ^ & f п/ a+1 I ab :а\\ ~

P(£ | y ) P(x I y )

T ZI/& \\ f P(x | X ) Х-1 N r a+1|ae о/ а+1 а+1 \\ i a+1 _

Таким образом, соотношение (26) ~f -——/ . =/ j °\\x -Zj )dx =

может быть представлено в в виде: P(x | y )

„_ N . P( | x" )

p^\\y ) = =1н j|a P(;~:+1 |yab:)■ (34)

Р(х"+ 1 | г" )Р(Х"+1 | уа :а )

= Р(Г | уа :а )| -——ь-¿ха+\\ (29) Подставив теперь выражения из соотР(х"+ \\уа"а) ношений (17), (33) и (34) в формулу (29), погде знаменатель подынтегрального выраже- лучаем приближенное равенство (32). Полания рассчитывается по формуле (12). гается, что Р(ха1 | ха) = Р((Га) 1 (ха+1,ха)), где

Рассмотрим маргинальное сглажи- (Га )-1 — обратная к Га функция, которая возвание с использованием частиц (Нигее1ег, вращает значение а по заданным Г+1 ,Г\\ см.

КипвеЬ, 1998; СоЖШ е1 а1., 2004; Багкка, (4). Описанный метод имеет вычислительную

через изменение весов частиц, поскольку сво- 1998; Godsi11 е1 а1., 2004; Багкка, 2013).

РЕЗУЛЬТАТЫ

Для выполнения расчетов по предложенной методике пакет Fishmetica в среде Julia был пополнен функциями для реализации фильтрации, сглаживания и предсказания временных рядов численности и уловов по возрастам с применением описанных в предыдущем разделе методов частиц. Оцененные численности каждой возраста характеризуется теперь не скалярными значением или параметрами известного распределения, а выборкой некоторого объема N из соответствующего распределения. По этим выборкам можно получать такие приближенные числовые характеристики, как средние, медианы, дисперсии и т.д., используя функции из статистического пакета StatBase.

Пробные расчеты проведены на первом из трех комплектов тестовых данных из ВНИРО для программ, применяемых в институтах отрасли при моделировании динамики запасов. Для этого комплекта в предположении о гауссовском распределении невязок результаты расчетов сравнивались со значениями параметров результирующих (также гауссовских) распределений для байесовских оценок, полученными с использованием явных соотношений из статьи И. И. Шевченко

Применение методики для отдельной когорты иллюстрируется на примере 36-й по порядку когорты, включающей рыб, которые в 1994 г. имели возраст один год (ab = 1, ae = 10). Вычисления выполнены в среде IJulia с использованием различных библиотек языка Julia, включая пакет функций и типов Fishmetica (https://github.com/ picestcode/Fishmetica.jl).

На рисунке 1 исходная динамика логарифмов численности изображена сплошной линией, отфильтрованные значения математических ожиданий численности и предсказания на три года, рассчитанные по формулам из статьи И. И. Шевченко (2019),— коротким пунктиром, а по аппроксимирующей выборке (при N = 10000) — длинным пунктиром. Соответствующие дисперсии и их оценки показаны на рисунке 2. Плотность распределения и гистограмма отфильтрованРис. 2. Фильтрация и предсказание: исходные и средние логарифмы численности для 36-й когорты по возрастам.

ных логарифмов численности для возраста с максимальной ошибкой даны на рисунке 3.

На рисунках 5—7 приведены аналогичные величины для сглаживания.

В таблицах 1—3 собраны характеристики вычислительных процессов для случая нормальных распределений при использовании различного количества частиц N для 36-й когорты. Для каждого N приведены приближенные значения времени ( t ), затраченного на расчеты параметров (сек) на виртуальной машине Intel (R) Xeon (R) CPU E5—2660 0 @ 2.20GHz с 32 Gb памяти в Центре обработки и хранения данных Тихоокеанского филиала ВНИРО (ФГБНУ «ТИНРО »), минимальное эффективное число частиц ( minne ), нормы разностей математического ожидания ( |Л|| ) и дисперсии ( Л2 ) от значений, полученных по формулам из статьи И. И. Шевченко (2019). Отметим, что получаемые результаты зависят от конкретных реализаций выборок, которые генерируются в процессе вычислений с применением генераторов случайных выборок. Полученные результаты показывают, что

увеличение числа частиц порой приводит к значительному росту времени расчетов, но необязательно приводит к увеличению точности расчетов одновременно всех параметров гауссовских распределений.

На рисунке 8 показаны суммарные уловы по годам для исходных данных (сплошная линия), а также рассчитанные по соотношениям из статьи И. И. Шевченко (2019) (короткий пунктир) и с использованием метода частиц при N = 5000 (длинный пунктир). Время расчетов по данным для всех 45 составило 70402.68 сек на упомянутой виртуальной машине с 32 процессорами и использованием библиотеки Distributed для распараллеливания вычислений.

ОБСУЖДЕНИЕ

Извлечение информации из времен -ных рядов численности и уловов по возрастам может проводиться с использованием моделей, построенных с применением уравнений когорт и уловов, а также дополнительных предположений о зависимостях, которые опиРис. 3. Фильтрация и предсказание: дисперсии логарифмов численности для численности 36-й когорты в возрасте 7 лет.

Рис. 4. Фильтрация: функция плотности вероятности и нормированная гистограмма для численности 36-й когорты в возрасте 7 лет.

Рис. 5. Предсказание: функция плотности вероятности и нормированная гистограмма для оценки численности 36-й когорты через 3 года.

Рис. 6. Сглаживание: исходные и средние логарифмы численности для 36-й когорты по возрастам.

Рис. 7. Сглаживание: дисперсии логарифмов численности для 36-й когорты по возрастам. ВОПРОСЫ РЫБОЛОВСТВА том 21 №2 2020 245

Рис. 8. Сглаживание: функция плотности вероятности и нормированная гистограмма для 36-й когорты в возрасте 7 лет.

Рис. 9. Сглаживание: суммарные уловы по годам.

Таблица 1. Характеристики вычислительных процессов для фильтрации

N t minne А|| Аа2

Таблица 2. Характеристики вычислительных процессов для предсказания

N t min ne 1 а

Таблица 3. Характеристики вычислительных процессов для сглаживания

N t min ne 1 а

сывают коэффициенты естественной и промысловой смертности, и о виде распределений невязок (Шевченко, 2017). Тогда логарифмы численности рассматриваются в качестве ненаблюдаемого состояния, а логарифмы уловов выступают в качестве наблюдений. Вместо точечных значений анализируются плотности распределения соответствующих случайных величин, которые позволяют вычислять любые числовые характеристики распределений (средние, медианы, дисперсии и т.д.). Байесовское оценивание различных апостериорных плотностей может быть сведено к последовательному решению некоторых уравнений. По этим уравнениям аналитические описания плотностей получаются только для очень специальных случаев (линейных гауссовских моделей) (Шевченко, 2019). В общем случае требуется применение численных методов.

В работе рассмотрено использование метода частиц (Монте-Карло) для проведения байесовского оценивания апостериорных плотностей. Этот метод применим даже в случаях нелинейных и/или негауссовских моделей динамики когорт. Пакет Fishmetica в среде Julia пополнен функциями, предназначенными для последовательной генерации выборок при проведении фильтрации, предсказания и сглаживания. Точность методов частиц при конечных объемах выборок достаточно сложно оценивать (Crisan, Doucet, 2019). В качестве тестового выбран первый набор данных ВНИРО для программ, применяемых в институтах отрасли при моделировании динамики запасов. Для расчетов по этому набору в предположении о нормальности невязок можно пользоваться также явными соотношениями из статьи

И. И. Шевченко (2019). Путем сравнения со значениями, вычисленными с использованием известных формул, продемонстрировано влияние количества частиц, используемых при генерации выборок, на точность оценки параметров результирующих апостериорных плотностей. Приведены также оценки минимального количества эффективных частиц на отдельных шагах вычислений. В целом увеличение числа частиц при одноразовых вычислениях в диапазоне от 1000 до 10000 не позволило гарантировано увеличить точность определения сразу обоих параметров известных гауссов-ских распределений для отдельной когорты (табл. 1—3). Хотя сглаженные оценки математических ожиданий отслеживают основные тенденции, которые наблюдаются в исходных данных (рис. 5, 8).

Продолжение исследований связано с анализом временных рядов численности и уловов по возрастам в предположении о лапласовском распределении невязок, при котором апостериорные распределения уже не являются лапласовскими. Поэтому лапла-совское распределение с технической точки зрения менее приспособлено для байесовского оценивания (Russell, Norvig, 2010). Однако оно позволяет в большей мере учесть особенности временных рядов реальных наблюдений по сравнению с гауссовским (Kotz et al., 2001).

СПИСОК ЛИТЕРАТУРЫ

Шевченко И. И. Моделирование промысловых запасов при известных оценках возрастной структуры популяций и уловов // Вопр. рыболовства. 2017. Т. 18. № 4. С. 507519.

Шевченко И. И. Моделирование промысловых запасов при известных оценках возрастной структуры популяций и уловов. II // Вопр. рыболовства. 2019. Т. 20. № 2. С. 152-163.

Aldrin M., Aanes S., Subbey S. Comments on incongruous formulations in the sam (state-space assessment model) model and consequences for fish stock assessment // Fisheries

Research, 2019. № 210. P. 224-227.

Crisan D., Doucet A. A survey of convergence results on particle filtering methods for practitioners / / IEEE Transactions on signal processing. 2002. V. 50. № 3. P. 736-746.

Doucet A., Godsill S., Andrieu C. On sequential monte carlo sampling methods for bayesian filtering // Statistics and computing.

Godsill S. J., Doucet A., West M. Monte carlo smoothing for nonlinear time series // J. American statistical association,

Gordon N., Salmond D. J.,

Smith A. F. M. Novel approach to nonlinear/ non-gaussian bayesian state estimation // IEE Proceedings F (Radar and Signal Processing). 1993. V. 140. P. 107-113.

Hurzeler M., Kunsch H. R. Monte carlo approximations for general state-space models // J. Computational and Graphical Statistics. 1998. V. 7. № 2. P. 175-193.

Kitagawa G. Monte carlo filter and smoother for non-gaussian nonlinear state space models // J. computational and graphical statistics. 1996. V. 5. № 1. P. 1-25.

Kotz S., Kozubowski T. J., Podgor-ski K. The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance. Basel: Birkhauser. 2001. 349 p.

Nielsen A, Berg C. W. Response to: Comments on incongruous formulations in the sam (state-space assessment model) model and consequences for fish stock assessment // Fisheries Research. 2019. № 210. P. 228229.

Russell S. J., Norvig P. Artificial Intelligence: A Modern Approach. Upper Saddle River: Prentice Hall, 2010. 1152 p.

Sarkka S. Bayesian filtering and smoothing. Cambridge: Cambridge University Press,

Schnute J. T. A general framework for developing sequential fisheries models // Canadian J. Fisheries and Aquatic Sciences. 1994. V. 51. № 8. P. 1676-1688.

Smith A. F. M., Gelfand A. E. Bayes- Valpine P. de, Hilborn R. State-space

ian statistics without tears: a sampling-resam- likelihoods for nonlinear fisheries time-series // pling perspective // The American Statistician. Canadian J. Fisheries and Aquatic Sciences.

SEQUENTIAL PARTICLE METHODS FOR BAYESIAN EVALUATIONS OF ABUNDANCE AND CATCH AT AGE

© 2020 y I. I. Shevchenko

The Pacific Branch Russian Federal Research Institute of Fisheries and Oceanography,

Vladivostok, 690091

A cohort population dynamics may be represented as a hidden Bayesian model with abundances as hidden states and catches as observations. Using these models, one can evaluate posterior densities and calculate such point-wise characteristics as means, medians, variances and so on. With rare exceptions (as linear Gaussian models), the recurrence equations met by the posterior densities have no analytic solutions. We describe several particle (Monte Carlo) methods that may be used for the density approximations and evaluations of their statistical quantities for nonlinear non-Gaussian models as well. The Fishmetica package was extended with functions for generating samples and masses for time series filtering, prediction, and smoothing. Evaluations in Julia were made for a test dataset with assumed Gaussian distributions of residuals and different numbers of particles. Numerical results were compared with known analytic ones.