Хотя я не знаком с внутренней работой алгоритма подгонки ZIP-HMM, есть некоторые очевидные различия в том, что вы реализовали в модели Стэна, и в том, как описывает себя алгоритм оптимизации ZIP-HMM. Рассмотрение этих вопросов кажется достаточным для получения аналогичных результатов.

Различия между моделями

Начальная вероятность состояния

Значения, которые оценивает ZIP-HMM, в частности fit1$prior, указывают на то, что он включает в себя возможность узнать вероятность для начального состояния. Однако в модели Стэна это фиксированное значение 1:1.

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

Это следует изменить, чтобы позволить модели оценить начальное состояние.

Приоры по параметрам (опционально)

Модель Стэна имеет неплоские априорные значения для lambda и theta, но, по-видимому, ZIP-HMM не взвешивает полученные значения. Если кто-то хотел более реалистично имитировать ZIP-HMM, то плоские приоры были бы лучше. Однако возможность иметь неплоские априорные значения в Стэне на самом деле дает возможность разработать более хорошо настроенную модель, чем это достижимо с помощью стандартных алгоритмов вывода HMM.

Нулевая инфляция в штате 1

Из документации метода fasthmmfit

Алгоритм быстрого градиентного спуска/стохастического градиентного спуска для изучения параметров в специализированной скрытой марковской модели с нулевым раздуванием, где нулевое раздувание происходит только в состоянии 1. [курсив добавлен]

Модель Стэна предполагает нулевую инфляцию во всех штатах. Вероятно, поэтому оценочное значение theta занижено по сравнению с оценкой ZIP-HMM MAP.

Государственный заказ

При оценке дискретных скрытых состояний или кластеров в Stan можно использовать упорядоченный вектор в качестве трюка, чтобы устранение проблем с переключением ярлыков. Здесь это эффективно достигается с помощью

positive_ordered[m] lambda;

Однако, поскольку ZIP-HMM имеет нулевую инфляцию только в первом состоянии, правильная реализация этого поведения в Стэне требует предварительного знания того, каков ранг lambda для «первого» состояния. Это кажется очень проблематичным для обобщения этого кода. А пока давайте просто двигаться вперед, предполагая, что мы всегда можем каким-то образом восстановить эту информацию. В этом конкретном случае мы предположим, что состояние 1 в HMM имеет более высокое значение lambda и, следовательно, будет состоянием 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, мы сначала компилируем ее и запускаем оценку 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

показывает, что они точно отражают значения, которые вывел 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