r - ZIP - 隐马尔可夫模型 r Stan

Question

我正在尝试与 Stan 一起调整零膨胀泊松隐马尔可夫模型。对于过去论坛中的 Poisson-HMM，显示了此设置。 see link .

同时调整ZIP与经典理论有据可查的代码和模型。

---

zip

library(ziphsmm)
set.seed(123)
prior_init <- c(0.5,0.5)
emit_init <- c(20,6)
zero_init <- c(0.5,0)
tpm <- matrix(c(0.9, 0.1, 0.2, 0.8),2,2,byrow=TRUE)
result <- hmmsim(n=100,M=2,prior=prior_init, tpm_parm=tpm,emit_parm=emit_init,zeroprop=zero_init)
y <- result$series
serie <- data.frame(y = result$series, m = result$state)

fit1 <-  fasthmmfit(y,x=NULL,ntimes=NULL,M=2,prior_init,tpm,
                    emit_init,0.5, hessian=FALSE,method="BFGS", 
                    control=list(trace=1))
fit1
$prior
            [,1]
[1,] 0.997497445
[2,] 0.002502555

$tpm
          [,1]       [,2]
[1,] 0.9264945 0.07350553
[2,] 0.3303533 0.66964673

$zeroprop
[1] 0.6342182

$emit
          [,1]
[1,] 20.384688
[2,]  7.365498

$working_parm
[1] -5.9879373 -2.5340475  0.7065877  0.5503559  3.0147840  1.9968067

$negloglik
[1] 208.823

斯坦

library(rstan)

ZIPHMM <- 'data {
    int<lower=0> N;
    int<lower=0> y[N];
    int<lower=1> m;
  }

parameters {
    real<lower=0, upper=1> theta; //
    positive_ordered[m] lambda; //
    simplex[m] Gamma[m]; // tpm
  }

model {
  vector[m] log_Gamma_tr[m];
  vector[m] lp;
  vector[m] lp_p1;

  // priors
  lambda ~ gamma(0.1,0.01);
  theta ~ beta(0.05, 0.05);

  // transposing tpm and taking the log of each entry
  for(i in 1:m)
    for(j in 1:m)
      log_Gamma_tr[j, i] = log(Gamma[i, j]);

  lp = rep_vector(-log(m), m); // 

    for(n in 1:N) {
      for(j in 1:m){
        if (y[n] == 0)
          lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp) +
                     log_sum_exp(bernoulli_lpmf(1 | theta),
                     bernoulli_lpmf(0 | theta) + poisson_lpmf(y[n] | lambda[j]));
        else
          lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp) + 
                     bernoulli_lpmf(0 | theta) + 
                     poisson_lpmf(y[n] | lambda[j]);
                   }
      lp = lp_p1;
                  }
    target += log_sum_exp(lp);
}'
mod_ZIP <- stan(model_code = ZIPHMM, data=list(N=length(y), y=y, m=2), iter=1000, chains=1)
print(mod_ZIP,digits_summary = 3)
               mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff  Rhat
theta         0.518   0.002 0.052    0.417    0.484    0.518    0.554    0.621   568 0.998
lambda[1]     7.620   0.039 0.787    6.190    7.038    7.619    8.194    9.132   404 1.005
lambda[2]    20.544   0.039 0.957   18.861   19.891   20.500   21.189   22.611   614 1.005
Gamma[1,1]    0.664   0.004 0.094    0.473    0.604    0.669    0.730    0.841   541 0.998
Gamma[1,2]    0.336   0.004 0.094    0.159    0.270    0.331    0.396    0.527   541 0.998
Gamma[2,1]    0.163   0.003 0.066    0.057    0.114    0.159    0.201    0.312   522 0.999
Gamma[2,2]    0.837   0.003 0.066    0.688    0.799    0.841    0.886    0.943   522 0.999
lp__       -222.870   0.133 1.683 -227.154 -223.760 -222.469 -221.691 -220.689   161 0.999

真值

real = list(tpm = tpm, 
     zeroprop = nrow(serie[serie$m == 1 & serie$y == 0, ]) / nrow(serie[serie$m == 1,]),
     emit = t(t(tapply(serie$y[serie$y != 0],serie$m[serie$y != 0], mean))))
real
$tpm
     [,1] [,2]
[1,]  0.9  0.1
[2,]  0.2  0.8

$zeroprop
[1] 0.6341463

$emit
       [,1]
1 20.433333
2  7.277778

估计给某人很奇怪可以帮助我知道我做错了。正如我们看到的，stan zeroprop 的估计值 = 0.518，而实际值为 0.634，另一方面，t.p.m. 的值。在 stan 中，它们相距很远，lambda1 = 7.62 和 lambda2 = 20.54 的平均值虽然它们足够接近，但与实际的 20.43 和 7.27 的顺序不同。我认为我在 Stan 中定义模型时犯了一些错误，但我不知道是哪个错误。

Answer 1

虽然我不知道 ZIP-HMM 拟合算法的内部工作原理，但您在 Stan 模型中实现的内容与 ZIP-HMM 优化算法如何描述自身存在一些明显差异。解决这些问题似乎足以产生类似的结果。

模型之间的差异

初始状态概率

ZIP-HMM 估计的值，特别是 fit1$prior，表明它具有学习初始状态概率的能力。然而，在 Stan 模型中，这是固定为 1:1

lp = rep_vector(-log(m), m);

这应该更改为允许模型估计初始状态。

参数先验(可选)

Stan 模型在 lambda 和 theta 上有不平坦的先验，但大概 ZIP-HMM 没有对它到达的特定值进行加权。如果想更真实地模仿 ZIP-HMM，那么平坦的先验会更好。然而，在 Stan 中拥有非平坦先验的能力确实是一个机会，可以开发比标准 HMM 推理算法更优化的模型。

状态 1 的零通胀

来自fasthmmfit方法的文档

Fast gradient descent / stochastic gradient descent algorithm to learn the parameters in a specialized zero-inflated hidden Markov model, where zero-inflation only happens in State 1. [emphasis added]

Stan 模型假设所有州的通货膨胀率为零。这可能是估计的 theta 值相对于 ZIP-HMM MAP 估计值缩小的原因。

状态排序

在估计 Stan 中的离散潜在状态或集群时，可以使用有序向量作为 mitigate against label switching issues 的技巧.这是通过

有效实现的

positive_ordered[m] lambda;

但是，由于 ZIP-HMM 仅在第一状态具有零膨胀，因此在 Stan 中正确实现此行为需要事先了解 lambda 的等级对于“第一”状态.这对于泛化这段代码来说似乎很有问题。现在，假设我们总能以某种方式恢复这些信息，让我们继续前进。在这种特定情况下，我们将假设 HMM 中的状态 1 具有更高的 lambda 值，因此将是 Stan 模型中的状态 2。

更新的 Stan 模型

在模型中加入上述变化应该是这样的

斯坦模型

data {
  int<lower=0> N;    // length of chain
  int<lower=0> y[N]; // emissions
  int<lower=1> m;    // num states
}

parameters {
  simplex[m] start_pos;         // initial pos probs
  real<lower=0, upper=1> theta; // zero-inflation parameter
  positive_ordered[m] lambda;   // emission poisson params
  simplex[m] Gamma[m];          // transition prob matrix
}

model {
  vector[m] log_Gamma_tr[m];
  vector[m] lp;
  vector[m] lp_p1;

  // transposing tpm and taking the log of each entry
  for (i in 1:m) {
    for (j in 1:m) { 
      log_Gamma_tr[j, i] = log(Gamma[i, j]);
    }
  }

  // initial position log-lik
  lp = log(start_pos);

  for (n in 1:N) {
    for (j in 1:m) {
      // log-lik for state
      lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp);

      // log-lik for emission
      if (j == 2) { // assuming only state 2 has zero-inflation
        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; // log-lik for next position
  }
  target += log_sum_exp(lp);
}

map 估计

将上面的内容作为字符串变量 code.ZIPHMM 加载，我们首先编译它并运行 MAP 估计(因为 MAP 估计的行为最像 HMM 拟合算法):

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

显示它们密切反射(reflect)了 fasthmmfit 推断的值，除了状态顺序被切换。

后验采样

该模型也可以与 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