混合 Pareto 和 Normal Stan 模型不起作用答案

【问题标题】：Mixed Pareto and Normal Stan model not working混合 Pareto 和 Normal Stan 模型不起作用
【发布时间】：2015-02-07 14:32:07
【问题描述】：

我正在尝试通过 rstan 学习 Stan（因为我熟悉 R）。我试过运行一个简单的混合 Pareto 和 Normal 模型。它编译得很好（据我所知），但它无法采样，给我错误：

"(-2, 2) 之间的初始化在 100 次尝试后失败。尝试指定初始值、缩小约束值范围或重新参数化模型。

调用采样器时出错；采样未完成”

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

我的 R+rstan 代码如下：

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 拟合 Pareto 模型的例子，所以我很难交叉验证我的方法。 p>

【问题讨论】：

仍然不采样的实验模型代码示例（超限制 +ve 先验接近正确解）：HMF.stan N; real x[N]; real sdev[N]; } 参数 { real y_min; real alpha; real xtrue[N]; } model { y_min ~ normal (10, 5) T[5,15]; alpha ~ normal(2, 5) T[1,3]; xtrue ~ normal(y_min, alpha); for(i in 1:N){ x[i] ~正常(xtrue[i], sdev[i]); } } "

标签： r jags stan rstan

【解决方案1】：

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

我能够使用模型在 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 打开一个错误，指出 x ~ pareto(y_min, alpha) 生成低于 y_min 的值。~~

代码适用于最新的 Stan 版本。请先升级，这个bug好像已经修复了。

【讨论】：

据我所知，我使用的是最新的 R (3.1.2) / rstan (2.6.0) / stan。不过，我正在使用 Mountain Lion (10.8.0)。你的设置是什么？
我刚刚从 git 获取了最新版本的 Stan 和 CmdStan。

【解决方案2】：

根本问题是参数上声明的支持不匹配 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 ,无穷大）

最简单的做法是将参数声明为 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]; 也许均匀分布或四参数 beta 分布比截断对数正态分布更合适。

最后，虽然逻辑上没有错 for(i in 1:N){ x[i] ~ normal(xtrue[i], sdev[i]); } 在计算上比逻辑上等价的语句差得多 x ~ normal(xtrue, sdev);

【讨论】：

如果我应用非常严格的范围：参数 { real y_min;真正的阿尔法；实 xtrue[N]; } 模型 { y_min ~ lognormal(1, 1) T[9,11]; alpha ~ lognormal(1, 1) T[1, 3]; xtrue ~ pareto(y_min, alpha); x ~ 正常（xtrue，sdev）；我打印 y_min, alpha 和 xtrue 我可以看到它在以下情况下出错：y_min=11 alpha=1 xtrue=[inf, 8.02369e+219, 1.23259e+48, inf, 2.94639e+23, inf, 11, 5.25037 e+21, 11, 11] 即，当它在数字限制处采样时，它看起来像失败了。其他发行版不会发生这种情况 - 正常情况下是可以的。

【解决方案3】：

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

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

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

【讨论】：