added HMM simultations and new section in manual in time-series chapter

Author: Bob Carpenter

misc/hmm/hmm-analytic.stan

@@ -0,0 +1,31 @@
+data {
+  int<lower=1> K;  // num categories
+  int<lower=1> V;  // num words
+  int<lower=0> T;  // num instances
+  int<lower=1,upper=V> w[T]; // words
+  int<lower=1,upper=K> z[T]; // categories
+  vector<lower=0>[K] alpha;  // transit prior
+  vector<lower=0>[V] beta;   // emit prior
+}
+transformed data {
+  vector<lower=0>[K] alpha_post[K];
+  vector<lower=0>[V] beta_post[K];
+  for (k in 1:K) 
+    alpha_post[k] <- alpha;
+  for (t in 2:T)
+    alpha_post[z[t-1],z[t]] <- alpha_post[z[t-1],z[t]] + 1;
+  for (k in 1:K)
+    beta_post[k] <- beta;
+  for (t in 1:T)
+    beta_post[z[t],w[t]] <- beta_post[z[t],w[t]] + 1;
+}
+parameters {
+  simplex[K] theta[K];  // transit probs
+  simplex[V] phi[K];    // emit probs
+}
+model {
+  for (k in 1:K) 
+    theta[k] ~ dirichlet(alpha_post[k]);
+  for (k in 1:K)
+    phi[k] ~ dirichlet(beta_post[k]);
+}

misc/hmm/hmm-fit-semisup.stan

@@ -0,0 +1,73 @@
+data {
+  int<lower=1> K;  // num categories
+  int<lower=1> V;  // num words
+  int<lower=0> T;  // num supervised items
+  int<lower=1> T_unsup;  // num unsupervised items
+  int<lower=1,upper=V> w[T]; // words
+  int<lower=1,upper=K> z[T]; // categories
+  int<lower=1,upper=V> u[T_unsup]; // unsup words
+  vector<lower=0>[K] alpha;  // transit prior
+  vector<lower=0>[V] beta;   // emit prior
+}
+parameters {
+  simplex[K] theta[K];  // transit probs
+  simplex[V] phi[K];    // emit probs
+}
+model {
+  for (k in 1:K) 
+    theta[k] ~ dirichlet(alpha);
+  for (k in 1:K)
+    phi[k] ~ dirichlet(beta);
+  for (t in 1:T)
+    w[t] ~ categorical(phi[z[t]]);
+  for (t in 2:T)
+    z[t] ~ categorical(theta[z[t-1]]);
+
+  { 
+    // forward algorithm computes log p(u|...)
+    real acc[K];
+    real gamma[T_unsup,K];
+    for (k in 1:K)
+      gamma[1,k] <- log(phi[k,u[1]]);
+    for (t in 2:T_unsup) {
+      for (k in 1:K) {
+        for (j in 1:K)
+          acc[j] <- gamma[t-1,j] + log(theta[j,k]) + log(phi[k,u[t]]);
+        gamma[t,k] <- log_sum_exp(acc);
+      }
+    }
+    lp__ <- lp__ + log_sum_exp(gamma[T_unsup]);
+  }
+}
+generated quantities {
+  int<lower=1,upper=K> y_star[T_unsup];
+  real log_p_y_star;
+  { 
+    // Viterbi algorithm
+    int back_ptr[T_unsup,K];
+    real best_logp[T_unsup,K];
+    real best_total_logp;
+    for (k in 1:K)
+      best_logp[1,K] <- log(phi[k,u[1]]);
+    for (t in 2:T_unsup) {
+      for (k in 1:K) {
+        best_logp[t,k] <- negative_infinity();
+        for (j in 1:K) {
+          real logp;
+          logp <- best_logp[t-1,j] + log(theta[j,k]) + log(phi[k,u[t]]);
+          if (logp > best_logp[t,k]) {
+            back_ptr[t,k] <- j;
+            best_logp[t,k] <- logp;
+          }
+        }
+      }
+    }
+    log_p_y_star <- max(best_logp[T_unsup]);
+    for (k in 1:K)
+      if (best_logp[T_unsup,k] == log_p_y_star)
+        y_star[T_unsup] <- k;
+    for (t in 1:(T_unsup - 1))
+      y_star[T_unsup - t] <- back_ptr[T_unsup - t + 1, 
+                                      y_star[T_unsup - t + 1]];
+  }
+}

misc/hmm/hmm-semisup.stan

@@ -0,0 +1,41 @@
+library("MCMCpack")
+library("rstan")
+
+# CONSTANTS
+K <- 3;
+V <- 10;
+T <- 100;
+T_unsup <- 500;
+alpha <- rep(1,K);
+beta <- rep(0.1,V);
+
+# DATA
+w <- rep(0,T);
+z <- rep(0,T);
+u <- rep(0,T_unsup);
+
+# PARAMETERS
+theta <- rdirichlet(K,alpha);
+phi <- rdirichlet(K,beta);
+
+# SIMULATE DATA
+
+# supervised
+z[1] <- sample(1:K,1);
+for (t in 2:T)
+  z[t] <- sample(1:K,1,replace=TRUE,theta[z[t - 1], 1:K]);
+for (t in 1:T)
+  w[t] <- sample(1:V,1,replace=TRUE,phi[z[t],1:V]);
+
+# unsupervised
+y <- rep(0,T_unsup);  
+y[1] <- sample(1:K,1);
+for (t in 2:T_unsup)
+  y[t] <- sample(1:K,1,replace=TRUE,theta[y[t-1],1:K]);
+for (t in 1:T_unsup)
+  u[t] <- sample(1:V,1,replace=TRUE,phi[y[t], 1:V]);
+
+
+fit <- stan('hmm-semisup.stan', # 'hmm-fit-semisup.stan',
+            data=list(K=K,V=V,T=T,T_unsup=T_unsup,M=M,w=w,z=z,u=u,alpha=alpha,beta=beta),
+            iter=200, chains=1, init=0);  # fit = fit // reuse model

misc/hmm/hmm-sim.R

@@ -0,0 +1,38 @@
+data {
+  int<lower=1> K;  // num categories
+  int<lower=1> V;  // num words
+  int<lower=0> T;  // num instances
+  int<lower=1,upper=V> w[T]; // words
+  int<lower=1,upper=K> z[T]; // categories
+  vector<lower=0>[K] alpha;  // transit prior
+  vector<lower=0>[V] beta;   // emit prior
+}
+transformed data {
+  int<lower=0> trans[K,K];
+  int<lower=0> emit[K,V];
+  for (k1 in 1:K) 
+    for (k2 in 1:K)
+      trans[k1,k2] <- 0;
+  for (t in 2:T)
+    trans[z[t - 1], z[t]] <- 1 + trans[z[t - 1], z[t]];
+  for (k in 1:K)
+    for (v in 1:V)
+      emit[k,v] <- 0;
+  for (t in 1:T)
+    emit[z[t], w[t]] <- 1 + emit[z[t], w[t]];
+}
+parameters {
+  simplex[K] theta[K];  // transit probs
+  simplex[V] phi[K];    // emit probs
+}
+model {
+  for (k in 1:K) 
+    theta[k] ~ dirichlet(alpha);
+  for (k in 1:K)
+    phi[k] ~ dirichlet(beta);
+
+  for (k in 1:K)
+    trans[k] ~ multinomial(theta[k]);
+  for (k in 1:K)
+    emit[k] ~ multinomial(phi[k]);
+}

misc/hmm/hmm-sufficient.stan

@@ -1,9 +1,9 @@
 data {
   int<lower=1> K;  // num categories
   int<lower=1> V;  // num words
-  int<lower=0> N;  // num instances
-  int<lower=1,upper=V> w[N]; // words
-  int<lower=1,upper=K> z[N]; // categories
+  int<lower=0> T;  // num instances
+  int<lower=1,upper=V> w[T]; // words
+  int<lower=1,upper=K> z[T]; // categories
   vector<lower=0>[K] alpha;  // transit prior
   vector<lower=0>[V] beta;   // emit prior
 }
@@ -16,8 +16,8 @@ model {
     theta[k] ~ dirichlet(alpha);
   for (k in 1:K)
     phi[k] ~ dirichlet(beta);
-  for (n in 1:N)
-    w[n] ~ categorical(phi[z[n]]);
-  for (n in 2:N)
-    z[n] ~ categorical(theta[z[n-1]]);
+  for (t in 1:T)
+    w[t] ~ categorical(phi[z[t]]);
+  for (t in 2:T)
+    z[t] ~ categorical(theta[z[t - 1]]);
 }