【问题标题】:stan programming for ETAS modelETAS 模型的 stan 编程
【发布时间】:2016-05-04 16:42:19
【问题描述】:

我是 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 代码:

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

我知道我没有将时间指定为向量。你能帮我在模型部分写下可能性吗?我面临写强度的问题。我认为我过去在 R 中写入时间 t[i] 的强度的方式不是在 STAN 中执行此操作的写入方式。

一小部分(仅包含20次)数据如下: DAT =列表(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) P>

【问题讨论】:

  • 如果要对参数使用统一先验,则必须声明其下限和上限。此时,您无需明确使用uniform() PDF。因此,您的参数块看起来像 real&lt;lower=0&gt; mu; real&lt;lower=1.1,upper=5&gt; p; real&lt;lower=1,upper=3&gt; c;

标签: bayesian stan rstan


【解决方案1】:

pow 函数目前不对向量或数组进行操作,因此您必须循环来构造强度。此外,我认为您的意思是将t 声明为长度为N 的实数数组,它看起来像real&lt;lower=0&gt; 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 中进行数值近似。

【讨论】:

  • 我必须对 (1:(j-1)) 中的每个 i 取 t[j] 和 t[i] 之间的差,然后将总和超过 1 / (t[j]-t [i]+c)^(p-1) 。我怎么能在 stan 中做到这一点?我在取款时遇到问题。
  • 写一个循环来计算总和。
  • @BobCarpenter 我编辑了我的问题。你能查一下吗?
  • Ben 向您展示了您需要一次进行一位求和的循环。将值 sto sum 放入临时数组并应用 sum() 函数会更有效,但在一切正常之前可能不值得这样做。
猜你喜欢
  • 2015-06-05
  • 2021-05-21
  • 2022-06-23
  • 2020-03-26
  • 1970-01-01
  • 1970-01-01
  • 2021-02-08
  • 2018-08-15
  • 1970-01-01
相关资源
最近更新 更多