statdis {RIPSeeker} | R Documentation |
Given a transition matrix A, returns the stationary distribution of a Markov chain by computing the eigen vectors of A.
statdis(A)
A |
Transition probability matrix, a squared matrix of probabilities (0 ≤ p ≤ 1) with row and column length equal to that of alpha and beta and row sum and column sum both equal to 1 (within some numerical deviation of 1e-6). |
w |
Stationary weights for the distributions of K components based on the transition probability matrix. |
Yue Li
Capp\'e, O. (2001). H2M : A set of MATLAB/OCTAVE functions for the EM estimation of mixtures and hidden Markov models. (http://perso.telecom-paristech.fr/cappe/h2m/)
# Simulate data TRANS_s <- matrix(c(0.9, 0.1, 0.3, 0.7), nrow=2, byrow=TRUE) alpha_s <- c(2, 4) beta_s <- c(1, 0.25) Total <- 100 x <- nbh_gen(TRANS_s, alpha_s, beta_s, Total); count <- x$count label <- x$label Total <- length(count) # dummy initialization TRANS0 <- matrix(rep(0.5,4), 2) alpha0 <- c(1, 20) beta0 <- c(1, 1) NIT_MAX <- 50 TOL <- 1e-100 nbh <- nbh_em(count, TRANS0, alpha0, beta0, NIT_MAX, TOL) map.accuracy <- length(which(max.col(nbh$postprob) == label))/Total vit <- nbh_vit(count, nbh$TRANS, nbh$alpha, nbh$beta) vit.accuracy <- length(which(vit$class == label))/Total # Plot the marginal distribution (in the stationnary regime) # Compute negative binomial distributions for all model states t <- 0:max(count) tmp <- nbh_em(t, nbh$TRANS, nbh$alpha, nbh$beta, 1) dens <- tmp$dens w <- statdis(nbh$TRANS) # Plot estimate of marginal probabilities marprob <- apply(t(dens) * (t(w) %*% matrix(1, ncol=length(t))), 2, sum) plot(t, marprob, pch=8, col="blue", main="Estimated marginal distribution") # Plot empirical estimated probabilities dhist <- matrix(0, ncol=length(t)) for(i in t){ dhist[1+i] <- sum(count == i)/Total } points(t, dhist, pch=3, col="red")