首页 文章

混合帕累托和普通斯坦模型不起作用

提问于
浏览 次
4

我正试图通过rstan学习斯坦（因为我熟悉R） . 我试过运行一个简单的混合Pareto和Normal模型 . 编译很好（据我所知），但它没有采样，给我错误：

“尝试100次后，（ - 2，2）之间的初始化失败 . 请尝试指定初始值，减少约束值的范围或重新参数化模型 .

调用采样器时发生错误;未完成抽样“

我只想说我已经尝试了各种方法来参数化，并尝试设置初始值，但都无济于事 .

我的Rrstan代码如下：

library(rstan)
rpareto = function(n, location, shape){location/runif(n)^(1/shape)}
sdvec=runif(1e3,0.1,1)
HMFtest=list(x=rpareto(1e3,10,2)+rnorm(1e3,0,sdvec), sdev=sdvec, N=1e3)

HMF.stan <- "
data {
  int<lower=0> N;
  real x[N];
  real sdev[N];
}
parameters {
  real<lower=0,upper=20> y_min;
  real<lower=0,upper=4> alpha;
  real xtrue[N];
}
model {
  y_min ~ lognormal(1, 1);
  alpha ~ lognormal(1, 1);
  xtrue ~ pareto(y_min, alpha);
  for(i in 1:N){
    x[i] ~ normal(xtrue[i], sdev[i]);
  }
}
"

stan.test <- stan(model_code=HMF.stan, data=HMFtest, pars=c('y_min','alpha'), chains=1, iter=30000, warmup=10000)

这个例子适用于JAGS（因此我也标记了JAGS）并且我可以发布代码是否有用 .

顺便说一句，如果我将帕累托分布更改为其他正态分布，它运行正常（当然，这给了我一个无意义的答案） .

关于我仍然认为JAGS不是Stan的任何建议，但是我找不到任何人与Stan拟合帕累托模型的例子，所以我很难交叉验证我的方法 .

r jags stan rstan

3 回答

  • 2

    错误消息意味着尝试的所有随机起始点产生零的可能性 .

    我能用模型在Stan中重现你的问题

    data {
  int<lower=0> N;
  real x[N];
  real sdev[N];
}
parameters {
  real<lower=0,upper=20> y_min;
  real<lower=0,upper=4> alpha;
  real xtrue[N];
}
model {
  y_min ~ lognormal(1, 1);
  alpha ~ lognormal(1, 1);
  print("y_min=", y_min, " alpha=", alpha);
  xtrue ~ pareto(y_min, alpha);
  print("xtrue: ", xtrue);
  x ~ normal(xtrue, sdev);
  print("x=", x);
}

    和数据

    N <- 6
sdev <- c(0.3339302,0.2936877,0.8540434,0.2399283,0.1014759,0.3717446)
x <- c(12.640112,10.502748,11.015629,29.382395,61.180509,12.772482)

    编译并运行Stan 2.0.1（现在已经很老了）我得到如下输出：

    y_min=4.49609:0 alpha=2.54906:0
xtrue: [0.992331:0,0.303142:0,0.180334:0,1.96009:0,0.903113:0,1.75711:0]
x=[12.6401,10.5027,11.0156,29.3824,61.1805,12.7725]
y_min=17.0143:0 alpha=1.67509:0
xtrue: [-1.40618:0,1.82026:0,1.67344:0,-0.973618:0,0.746502:0,1.93469:0]
x=[12.6401,10.5027,11.0156,29.3824,61.1805,12.7725]

    因此，虽然为y_min和alpha选择了合理的参数，但是帕累托生成的值也低于y_min . 在手册中，概率分布函数也不包含截断 . 我认为这是问题（用正态分布代替帕累托运行正常） . 我建议在github上用Stan打开一个bug，说明x~pareto（y_min，alpha）生成的值低于y_min .

    代码适用于最新的Stan版本 . 请先升级，这个bug似乎已经修复了一段时间 .

    回复于
  • 0

    基本问题是在参数 parameters { real<lower=0,upper=20> y_min; real<lower=0,upper=4> alpha; real xtrue[N]; } 上声明的支持与先验的样本空间之间的不匹配 model { y_min ~ lognormal(1, 1); alpha ~ lognormal(1, 1); xtrue ~ pareto(y_min, alpha); ...

    • y_min 被约束到（0,20）区间，但对数正态先验在整个正实线上传播质量单位

    • alpha 被约束到（0,20）区间，但是对数正态先验在整个正实线上传播质量单位

    • 最糟糕的是， xtrue 的每个元素都是无约束的 - 意味着它可以是整个实线上的任何东西---但是帕累托事先在区间内传播一个质量单位（y_min，Infinity）

    最简单的方法是将参数声明为 parameters { real<lower=0> y_min; real<lower=0> alpha; real<lower=y_min> xtrue[N]; } 原则上，您可以保持 y_min 和 alpha 的上限，并指定一些在声明的支持上集成到1的先验 . 一种粗略的方法是截断（将对数正态PDF除以未截断质量的数量）对数正态先验，如 model { y_min ~ lognormal(1, 1) T[,20]; alpha ~ lognormal(1, 1) T[,4]; 也许统一或四参数β分布比截断对数正态更合适 .

    最后，虽然它在逻辑上不是错误的 for(i in 1:N){ x[i] ~ normal(xtrue[i], sdev[i]); } 在计算上比逻辑等价的语句更糟糕 x ~ normal(xtrue, sdev);

    回复于
  • 2

    所有固定的人 . 事实证明，我的R，stan和c编译器的混合存在一些深层次的问题 .

    不想手工修理东西，我采取了大锤的方法，只是升级到Yosemite，然后从头开始安装关键组件 . 这似乎可以解决所有问题，现在我的链条收敛得非常好 .

    这是一个奇怪的，因为编译器可以构建rstan并编译/采样许多rstan示例 . 我不知道帕累托为什么会引起这些问题，但现在肯定已经解决了 .

    回复于

相关问题