虽然我不知道ZIP-HMM拟合算法的内部工作原理,但在Stan模型中实现的内容以及ZIP-HMM优化算法如何描述自身方面,还是存在一些明显的差异。解决这些问题似乎足以产生类似的结果。
模型之间的差异
初始状态概率
fit1$prior
,表示它包括学习初始状态概率的能力。然而,在斯坦模型中,这是固定的1:1
lp = rep_vector(-log(m), m);
这应该改变,以允许模型估计初始状态。
参数优先级(可选)
斯坦模型有非平坦的优先权
lambda
和
theta
,但据推测,ZIP-HMM没有对它到达的特定值进行加权。如果你想更真实地模仿ZIP-HMM,那么平面优先会更好。然而,与标准HMM推理算法相比,Stan中具有非平坦先验的能力确实是一个开发更优化模型的机会。
状态1零通货膨胀
fasthmmfit
方法
零通货膨胀只发生在状态1
. [着重强调]
斯坦模型假设所有州的通货膨胀率为零。这很可能就是为什么
θ
值相对于ZIP-HMM映射估计值进行压缩。
国家秩序
当估计Stan中离散的潜在状态或簇时,可以使用
有序向量
作为一个诡计
mitigate against label switching issues
. 这在这里可以有效地实现
positive_ordered[m] lambda;
然而,由于ZIP-HMM在第一个状态下只有零膨胀,所以在Stan中正确实现此行为需要事先知道
兰姆达
代表“第一”州。这对于泛化这段代码似乎很有问题。现在,让我们继续假设我们总能以某种方式恢复这些信息。在这种特定情况下,我们将假设HMM中的状态1具有更高的
兰姆达
值,因此在Stan模型中将是状态2。
更新的斯坦模型
在模型中合并上述更改应该类似于
斯坦模型
data {
int<lower=0> N;
int<lower=0> y[N];
int<lower=1> m;
}
parameters {
simplex[m] start_pos;
real<lower=0, upper=1> theta;
positive_ordered[m] lambda;
simplex[m] Gamma[m];
}
model {
vector[m] log_Gamma_tr[m];
vector[m] lp;
vector[m] lp_p1;
for (i in 1:m) {
for (j in 1:m) {
log_Gamma_tr[j, i] = log(Gamma[i, j]);
}
}
lp = log(start_pos);
for (n in 1:N) {
for (j in 1:m) {
lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp);
if (j == 2) {
if (y[n] == 0) {
lp_p1[j] += log_mix(theta, 0, poisson_lpmf(0 | lambda[j]));
} else {
lp_p1[j] += log1m(theta) + poisson_lpmf(y[n] | lambda[j]);
}
} else {
lp_p1[j] += poisson_lpmf(y[n] | lambda[j]);
}
}
lp = lp_p1;
}
target += log_sum_exp(lp);
}
地图估算
将上述内容作为字符串变量加载
code.ZIPHMM
model.ZIPHMM <- stan_model(model_code=code.ZIPHMM)
// note the use of some initialization on the params,
// otherwise it can occasionally converge to strange extrema
map.ZIPHMM <- optimizing(model.ZIPHMM, algorithm="BFGS",
data=list(N=length(y), y=y, m=2),
init=list(theta=0.5, lambda=c(5,10)))
检查估计参数
> map.ZIPHMM$par
start_pos[1] start_pos[2]
9.872279e-07 9.999990e-01
theta
6.342449e-01
lambda[1] lambda[2]
7.370525e+00 2.038363e+01
Gamma[1,1] Gamma[2,1] Gamma[1,2] Gamma[2,2]
6.700871e-01 7.253215e-02 3.299129e-01 9.274678e-01
显示它们紧密地反映了
快!快
推断,除了状态命令被切换。
从后面取样
这个模型也可以用MCMC来推断完全后验,
samples.ZIPHMM <- stan(model_code = code.ZIPHMM,
data=list(N=length(y), y=y, m=2),
iter=2000, chains=4)
> samples.ZIPHMM
Inference for Stan model: b29a2b7e93b53c78767aa4b0c11b62a0.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
start_pos[1] 0.45 0.00 0.29 0.02 0.20 0.43 0.69 0.97 6072 1
start_pos[2] 0.55 0.00 0.29 0.03 0.31 0.57 0.80 0.98 6072 1
theta 0.63 0.00 0.05 0.53 0.60 0.63 0.67 0.73 5710 1
lambda[1] 7.53 0.01 0.72 6.23 7.02 7.49 8.00 9.08 4036 1
lambda[2] 20.47 0.01 0.87 18.83 19.87 20.45 21.03 22.24 5964 1
Gamma[1,1] 0.65 0.00 0.11 0.43 0.57 0.65 0.72 0.84 5664 1
Gamma[1,2] 0.35 0.00 0.11 0.16 0.28 0.35 0.43 0.57 5664 1
Gamma[2,1] 0.08 0.00 0.03 0.03 0.06 0.08 0.10 0.16 5605 1
Gamma[2,2] 0.92 0.00 0.03 0.84 0.90 0.92 0.94 0.97 5605 1
lp__ -214.76 0.04 1.83 -219.21 -215.70 -214.43 -213.43 -212.25 1863 1