statdis: Returns the stationary distribution of a Markov chain.

Description Usage Arguments Value Author(s) References See Also Examples

Description

Given a transition matrix A, returns the stationary distribution of a Markov chain by computing the eigen vectors of A.

Usage

1

Arguments

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).

Value

w

Stationary weights for the distributions of K components based on the transition probability matrix.

Author(s)

Yue Li

References

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/)

See Also

nbh_em

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
	
# 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")

yueli-compbio/RIPSeeker documentation built on May 8, 2019, 2:34 a.m.