стандартное программирование для модели ETAS

Я новичок в STAN. Я работаю над временной моделью ETAS, моделью, используемой для моделирования землетрясений. Интенсивность землетрясения во время t [i] моделируется как:

h(t[i]|p,c,mu)=mu+sum((p-1)*(c^(p-1))*(1/((t[i]-t[1:(i-1)]+c)^(p-1))));

где t - время, а p, c, mu - три параметра. Я использую Rstan. Я написал для модели следующий стандартный код:

stan_etas="
data{
  int<lower=0> N;
  real<lower=0> t;
}
parameters{
  real<lower=0> mu;
  real<lower=1.005> p;
  real<lower=0> c;
}

Я знаю, что я не указывал время как вектор. Можете помочь написать правдоподобие в разделе модели? У меня проблема с записью интенсивности. Я думаю, что способ, которым я писал интенсивность в момент времени t [i] в ​​R, не является способом записи для этого в STAN.

Небольшая часть (содержащая только 20 раз) данных выглядит следующим образом: dat = list (0,0000,310,1907,948,4677,1007.2617,1029.7996,1065.7343,1199.8650, 1234.6809,1298.0234,1316.0350,1381.8400,1413.4311,1546.2059,1591.1326, 1669.5084, 1669.5084, 1738.9363,1745.5503,1797.9980,1895.6705,1936.3146)


rstan bayesian stan
person Geotas    schedule 04.05.2016    source источник
comment
Если вы собираетесь использовать предварительную униформу для параметра, вы должны объявить его нижнюю и верхнюю границы. И на этом этапе вам не нужно явно использовать uniform() PDF. Итак, ваш блок параметров будет выглядеть как real<lower=0> mu; real<lower=1.1,upper=5> p; real<lower=1,upper=3> c;   -  person Ben Goodrich    schedule 04.05.2016

Ответы (1)


arrow_upward
2
arrow_downward

Функция pow в настоящее время не работает с векторами или массивами, поэтому вам нужно выполнить цикл для построения интенсивности. Вдобавок, я думаю, вы хотели объявить t как реальный массив длины N, который будет выглядеть как real<lower=0> t[N];. Тогда в блоке модели у вас будет что-то вроде:

y[1] <- pow(c, -(p-1));
for (j in 2:N) {
  y[j] <- mu;
  for (i in 1:(j-1))
    y[j] <- y[j] + (p - 1) * c^(p-1) * 
            1 / (t[j]-t[i]+c)^(p-1);
 }

Однако в конечном итоге вам придется использовать функцию increment_log_prob() для регистрации логарифмической вероятности. Хотя я не знаком с моделью ETAS, документация пакета ETAS R утверждает, что он включает в себя интеграл, который в настоящее время не может быть вычислен численно в Stan.

person Ben Goodrich    schedule 04.05.2016
comment
Мне нужно взять разницу между t [j] и t [i] для каждого i в (1: (j-1)), а затем взять сумму по 1 / (t [j] -t [i] + c) ^ ( п-1). Как я мог это сделать? У меня проблема с получением суммы. - person Geotas; 05.05.2016
comment
Напишите цикл для вычисления суммы. - person Bob Carpenter; 05.05.2016
comment
@BobCarpenter Я отредактировал свой вопрос. Вы не могли бы проверить? - person Geotas; 05.05.2016
comment
Бен показывает вам цикл, в котором вам нужно посчитать сумму по частям. Немного эффективнее поместить значение sto sum во временный массив и применить функцию sum (), но, вероятно, не стоит этого делать, пока вы не все заработаете. - person Bob Carpenter; 07.05.2016