1
$\begingroup$

I'm trying to adjust a Zero Inflated Poisson Hidden Markov Model with Stan. For the Poisson-HMM in a past forum this setting was shown. see link.

While to adjust the ZIP with the classical theory is well documented the code and model.

Zero Inflation

It uses a parameter theta here there is a probability $\theta$ of drawing a zero, and a probability $1−\theta$ of drawing from Poisson($\lambda$)

$p(y_n | \theta, \lambda) = \begin{cases} \theta + (1-\theta)*e^{-\lambda} & \text{ if } y_n=0 \\ (1-\theta)*e^{-\lambda} \frac{\lambda^{y_n}}{y_n!} & \text{ if } y_n>0 \end{cases}$

Zero Inflation - Hidden Markov Model

Let $Y_t$ be the observed series and $C_t$ a homogeneous, unobserved Markov chain, withstate-space $E={e_1,· · ·, e_m}$. We suppose that, conditionally to $C_t=e_i$, $Y_t$ is distributed according to a ZIP of parameters $(\theta_i, \lambda_i)$.

$p(y_t | C_t ; \theta, \lambda) = \begin{cases} \theta_i + (1-\theta_i)*e^{-\lambda_i} & \text{ if } y_t=0 \\ (1-\theta_i)*e^{-\lambda_i} \frac{\lambda_i^{y_t}}{y_t!} & \text{ if } y_t>0 \end{cases}$

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
Stan 
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
True values 
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

Estimates give quite oddly to someone could help me to know that I am doing wrong.

I will expand my question, the proportion of zeros gives quite a distance with rStan theta = 0.518 and the real is 0.634, the same for the values of the transition matrix. Also the average of the values lambda1 = 7.62 and lambda2 = 20.54, while the real ones are lambda1 = 20.43 and lambda2 = 7.27. That is, they are crossed. I would expect to obtain estimates with rstan, similar to those of the ziphsmm package. I think I'm making some mistake in defining the model in Stan but I do not know which.

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
0
$\begingroup$

Making the appropriate adjustments for a ZIP-HMM of two states, where zeri inflation occurs only in state 1, we have the following code.

library(rstan)

code.ZIPHMM <- '
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_or;   // emission poisson params
  simplex[m] Gamma[m];          // transition prob matrix
}

transformed parameters{
  vector[m] lambda;
  lambda[1] = lambda_or[2];  
  lambda[2] = lambda_or[1];
}

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 == 1) { // assuming only state 1 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);
}'

samples.ZIPHMM <- stan(model_code = code.ZIPHMM, 
                       data=list(N=length(y), y=y, m=2), 
                       iter=1000, chains=1, 
                       pars = c("theta","lambda","Gamma","lp__"))

print(samples.ZIPHMM, digits_summary = 3)
               mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff  Rhat
theta         0.637   0.003 0.058    0.524    0.596    0.640    0.674    0.746   505 1.001
lambda[1]    20.422   0.033 0.869   18.824   19.794   20.411   20.985   22.102   696 0.999
lambda[2]     7.562   0.028 0.654    6.279    7.140    7.554    7.997    8.958   557 1.001
Gamma[1,1]    0.915   0.001 0.033    0.840    0.898    0.919    0.937    0.967   589 0.998
Gamma[1,2]    0.085   0.001 0.033    0.033    0.063    0.081    0.102    0.160   589 0.998
Gamma[2,1]    0.348   0.005 0.116    0.149    0.266    0.333    0.430    0.583   551 1.003
Gamma[2,2]    0.652   0.005 0.116    0.417    0.570    0.667    0.734    0.851   551 1.003
lp__       -214.914   0.104 1.866 -219.042 -216.092 -214.591 -213.451 -212.236   319 1.002

Now the estimates of the values with Stan match those of the ziphsmm package.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.