代码之家  ›  专栏  ›  技术社区  ›  Rafael Díaz

压缩隐马尔可夫模型

  •  -1
  • Rafael Díaz  · 技术社区  · 6 年前

    我试图用斯坦来调整一个零膨胀的Poisson隐马尔可夫模型。对于过去论坛中的Poisson HMM,显示了此设置。 see link .

    ZIP 用经典的理论很好地记录了代码和模型。

    enter image description here


    ziphsmm公司
    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
    

    给某人一个很奇怪的估计可以帮助我知道我做错了。如我们所见,当实际值为0.634时,stan zeroprop=0.518的估计值,另一方面,stan中的t.p.m.值非常遥远,平均值lambda1=7.62和lambda2=20.54,尽管它们以不同的顺序近似于实数20.43和7.27。我想我在定义Stan中的模型时犯了一些错误,但我不知道是哪个。

    0 回复  |  直到 6 年前
        1
  •  4
  •   merv    6 年前

    虽然我不知道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;    // 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);
    }
    

    地图估算

    将上述内容作为字符串变量加载 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