0

I am new in STAN. I am working on temporal ETAS model, a model used for modelling earthquakes.The intensity at earthquake occurrence time t[i] is modelled as-

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

where t is the time and p,c,mu are the three parameters. I am using the Rstan. I wrote the following stan code for the model:

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;
}

I know that I did not specified the time as vector. Can you help me to write the likelihood in the model section? I am facing problem writing the intensity. I think the way I used to write the intensity at time t[i] in R is not the write way to do this in STAN.

A small part (containing 20 times only) of the data is as follows: 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,1738.9363,1745.5503,1797.9980,1895.6705,1936.3146)

Geotas
  • 35
  • 6
  • If you are going to use a uniform prior on a parameter, then you must declare its lower and upper bounds. And at that point, you need not explicitly use the `uniform()` PDF. So, your parameters block would look like `real mu; real p; real c;` – Ben Goodrich May 04 '16 at 18:27

1 Answers1

2

The pow function currently does not operate over vectors or arrays, so you have to loop to construct the intensity. In addition, I think you meant to declare t as a real array of length N, which would look like real<lower=0> t[N];. Then in the model block, you would have something like:

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);
 }

However, you ultimately have to use the increment_log_prob() function to register the log-likelihood. Although I am not familiar with the ETAS model, the documentation of the ETAS R package claims it involves an integral, which currently cannot be approximated numerically in Stan.

Ben Goodrich
  • 4,870
  • 1
  • 19
  • 18
  • I have to take difference between t[j] and t[i] for each i in (1:(j-1)) and then take sum over 1 / (t[j]-t[i]+c)^(p-1) . How could I do that in stan? I am facing problem in taking the sum. – Geotas May 04 '16 at 22:54
  • Write a loop to compute the sum. – Bob Carpenter May 05 '16 at 02:26
  • @BobCarpenter I edited my question. Could you please check? – Geotas May 05 '16 at 08:37
  • Ben's showing you the loop you need to do the sum one bit at a time. It's a bit more efficient to put the value sto sum into a temp array and apply the sum() function, but probably not worth doing until you get everything working. – Bob Carpenter May 07 '16 at 01:33