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

Question

Я новичок в STAN. Я работаю на временной модели ETAS, модель, используемую для моделирования earthquakes.The интенсивности при землетрясении времени возникновения т [I] моделируется как-

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

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

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; 
} 
model{ 
    real y[N]; 
    //priors 
    mu~exponential(0.5); 
    p~uniform(1.1,5); 
    c~uniform(1,3); 
    //likelihood 
    y[1]<-pow(c,(-(p-1))); 
    for(j in 2:N){ 
    y[j]<-sum(pow(((t[j]-t[1:(j-1)])+c),(-(p-1)))); 
    } 
} 
" 
fit = stan(model_code=stan_etas, data=dat, iter=12000, warmup=2000, thin=10, chains=3) 

. Но я получаю это сообщение об ошибке следующим образом:

SYNTAX ERROR, MESSAGE(S) FROM PARSER: 

Indexed expression must have at least as many dimensions as number of  indexes supplied: 
indexed expression dimensionality = 0; indexes supplied = 1 

ERROR at line 20 

18: y[1]<-pow(c,(-(p-1))); 
19: for(j in 2:N){ 
20: y[j]<-sum(pow(((t[j]-t[1:(j-1)])+c),(-(p-1)))); 
         ^
21: } 

Error in stanc(file = file, model_code = model_code, model_name =  model_name, : 
    failed to parse Stan model '0bcc61f1791b1f924019273f75206514' due to the above error. 

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

Небольшая часть (содержащая 20 раз только) данных выглядит следующим образом: Дат = список (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,1738.9363,1745.5503,1797.9980,1895.6705,1936.3146)

РЕДАКТИРОВАТЬ: полный журнал правдоподобия код модели ETAS, как:
LL = сумма (журнал (мю + сумма ((р-1) (c^(p-1)) (1/((t [i] -t [1: (i-1)] + c)^(p-1)))))) - (mu * (p-1))))

; где tmax - максимальное время появления. Так что я редактировал код стан следующим образом: Дат = список (N = длина (ц), т = ц, Tmax = Tmax)

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

real h[N]; 
real a; 
real b[N]; 
real final; 
//priors 
mu~exponential(0.5); 
p~uniform(1.1,5); 
c~uniform(1,3); 
//likelihood 
for(j in 2:N){ 
    vector[j-1] y; 
    y[1]<-mu+pow(c,-(p-1)); 
    for(i in 1:(j-1)){ 
    y[i]<-((p - 1) * c^(p-1) * 
     1/(t[j]-t[i]+c)^(p-1)); 
    } 
    h[j]<-log(mu+sum(y)); 
} 
for (k in 1:N){ 
    b[k]<-1-((c^(p-1))/((tmax-t[k]+c)^(p-1))); 
} 
a<-mu*tmax; 
final<-sum(h)-a-sum(b); 
} 
" 
fit = stan(model_code=stan_etas, data=dat, iter=12000, warmup=2000, thin=10, chains=3) 

. Но я получаю эту ошибку:

"Rejecting initial value:"                    
[299] " Gradient evaluated at the initial value is not finite."            
[300] " Stan can't start sampling from this initial value."             
[301] "Initialization between (--2, -2) failed after 100 attempts. "           
[302] " Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model." 

Answer 1

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 package утверждает, что она включает в себя интеграл, который в настоящее время не может быть аппроксимирован численно в Stan.