-
Notifications
You must be signed in to change notification settings - Fork 43
/
Copy pathmain.R
143 lines (115 loc) · 4.45 KB
/
main.R
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
library(rstan)
library(shinystan)
source('common/R/math.R')
source('common/R/plots.R')
source('iohmm-reg/R/iohmm-sim.R')
# Set up ------------------------------------------------------------------
# Data
T.length = 300
K = 3
M = 4
R = 1
u.intercept = FALSE
w = matrix(
c(1.2, 0.5, 0.3, 0.1, 0.5, 1.2, 0.3, 0.1, 0.5, 0.1, 1.2, 0.1),
nrow = K, ncol = M, byrow = TRUE)
b = matrix(
c(5.0, 6.0, 7.0, 0.5, 1.0, 5.0, 0.1, -0.5, 0.1, -1.0, -5.0, 0.2),
nrow = K, ncol = M, byrow = TRUE)
s = c(0.2, 1.0, 2.5)
p1 = c(0.4, 0.2, 0.4)
# Markov Chain Monte Carlo
n.iter = 400
n.warmup = 200
n.chains = 1
n.cores = 1
n.thin = 1
n.seed = 9000
# Data simulation ---------------------------------------------------------
set.seed(n.seed)
u <- matrix(rnorm(T.length*M, 0, 1), nrow = T.length, ncol = M)
if (u.intercept)
u[, 1] = 1
dataset <- iohmm_sim(T.length, K, u, w, p1, obsmodel_reg, obs.pars = list(b = b, s = s))
# Data exploration --------------------------------------------------------
plot_inputoutput(x = dataset$x, u = dataset$u, z = dataset$z)
plot_inputprob(u = dataset$u, p.mat = dataset$p.mat, z = dataset$z)
# Model estimation --------------------------------------------------------
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
stan.model = 'iohmm-reg/stan/iohmm-reg.stan'
stan.data = list(
T = T.length,
K = K,
M = M,
u_tm = as.array(u),
x_t = dataset$x
)
stan.fit <- stan(file = stan.model,
data = stan.data, verbose = T,
iter = n.iter, warmup = n.warmup,
thin = n.thin, chains = n.chains,
cores = n.cores, seed = n.seed)
n.samples = (n.iter - n.warmup) * n.chains
# MCMC Diagnostics --------------------------------------------------------
options(digits = 2)
summary(stan.fit,
pars = c('p_1k', 'w_km', 'b_km', 's_k'),
probs = c(0.50))$summary
launch_shinystan(stan.fit)
# Extraction --------------------------------------------------------------
alpha_tk <- extract(stan.fit, pars = 'alpha_tk')[[1]]
gamma_tk <- extract(stan.fit, pars = 'gamma_tk')[[1]]
zstar_t <- extract(stan.fit, pars = 'zstar_t')[[1]]
hatx_t <- extract(stan.fit, pars = 'hatx_t')[[1]]
# Relabelling (ugly hack edition) -----------------------------------------
dataset$zrelab <- rep(0, T)
hard <- sapply(1:T.length, function(t, med) {
which.max(med[t, ])
}, med = apply(alpha_tk, c(2, 3),
function(x) {
quantile(x, c(0.50)) }))
tab <- table(hard = hard, original = dataset$z)
for (k in 1:K) {
dataset$zrelab[which(dataset$z == k)] <- which.max(tab[, k])
}
print("Label re-imputation (relabeling due to switching labels)")
table(new = dataset$zrelab, original = dataset$z)
# Estimation summary ------------------------------------------------------
print("Estimated initial state probabilities")
summary(stan.fit,
pars = c('p_1k'),
probs = c(0.10, 0.50, 0.90))$summary[, c(1, 3, 4, 5, 6)]
print("Estimated probabilities in the transition matrix")
summary(stan.fit,
pars = c('A_ij'),
probs = c(0.10, 0.50, 0.90))$summary[, c(1, 3, 4, 5, 6)]
print("Estimated regressors of hidden states")
summary(stan.fit,
pars = c('w_km'),
probs = c(0.10, 0.50, 0.90))$summary[, c(1, 3, 4, 5, 6)]
print("Estimated regressors and standard deviation of observations in each state")
summary(stan.fit,
pars = c('b_km', 's_k'),
probs = c(0.10, 0.50, 0.90))$summary[, c(1, 3, 4, 5, 6)]
print("Observations with no imputation by the smoother (check)")
sum(apply(gamma_tk, 1, rowSums) == 0)
# Inference summary -------------------------------------------------------
# Filtered and smoothed state probability plot
plot_stateprobability(alpha_tk, gamma_tk, 0.8, dataset$zrelab)
# Confusion matrix for hard (naive) classification
print("Estimated hidden states (hard naive classification using filtered probability)")
print(table(
estimated = apply(round(apply(alpha_tk, c(2, 3),
function(x) {
quantile(x, c(0.50)) })), 1, which.max),
real = dataset$zrelab))
# Jointly most likely state path (Viterbi decoding)
plot_statepath(zstar_t, dataset$zrelab)
# Confusion matrix for jointly most likely state path
print("Estimated hidden states for the jointly most likely path (Viterbi decoding)")
round(table(
actual = rep(dataset$zrelab, each = n.samples),
fit = zstar_t) / n.samples, 0)
# Fitted output
plot_outputfit(dataset$x, hatx_t, z = dataset$zrelab, TRUE)