BTS

hmm.r

library(HMM)

outputS = c('A', 'B', 'C', 'D', 'E')
hiddenS = c('S01', 'S02', 'S13', 'S14', 'S15', 'S56', 'S26', 'S37', 'S48', 'S68', 'E')
startP = c(0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10)
transP = rbind(
  t(c(0.00, 0.20, 0.20, 0.20, 0.20, 0.00, 0.00, 0.00, 0.00, 0.00, 0.20)),
  t(c(0.33, 0.00, 0.00, 0.00, 0.00, 0.00, 0.33, 0.00, 0.00, 0.00, 0.34)),
  t(c(0.20, 0.00, 0.00, 0.20, 0.20, 0.00, 0.00, 0.20, 0.00, 0.00, 0.20)),
  t(c(0.20, 0.00, 0.20, 0.00, 0.20, 0.00, 0.00, 0.00, 0.20, 0.00, 0.20)),
  t(c(0.20, 0.00, 0.20, 0.20, 0.00, 0.20, 0.00, 0.00, 0.00, 0.00, 0.20)),
  t(c(0.00, 0.00, 0.00, 0.00, 0.25, 0.00, 0.25, 0.00, 0.00, 0.25, 0.25)),
  t(c(0.00, 0.25, 0.00, 0.00, 0.00, 0.25, 0.00, 0.00, 0.00, 0.25, 0.25)),
  t(c(0.00, 0.00, 0.50, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50)),
  t(c(0.00, 0.00, 0.00, 0.33, 0.00, 0.00, 0.00, 0.00, 0.00, 0.33, 0.34)),
  t(c(0.00, 0.00, 0.00, 0.00, 0.00, 0.25, 0.25, 0.00, 0.25, 0.00, 0.25)),
  t(c(0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00))
)
emissP = rbind(
  t(c(1.00, 0.00, 0.00, 0.00, 0.00)),
  t(c(1.00, 0.00, 0.00, 0.00, 0.00)),
  t(c(0.00, 1.00, 0.00, 0.00, 0.00)),
  t(c(0.00, 1.00, 0.00, 0.00, 0.00)),
  t(c(0.00, 1.00, 0.00, 0.00, 0.00)),
  t(c(0.00, 0.00, 1.00, 0.00, 0.00)),
  t(c(0.00, 0.00, 1.00, 0.00, 0.00)),
  t(c(0.00, 0.00, 0.00, 1.00, 0.00)),
  t(c(0.00, 0.00, 0.00, 1.00, 0.00)),
  t(c(0.00, 0.00, 0.00, 1.00, 0.00)),
  t(c(0.00, 0.00, 0.00, 0.00, 1.00))
)

hmm <- initHMM(hiddenS, outputS, startP, transP, emissP)

obs_1 <- c('D', 'C', 'B', 'A', 'E')
obs_2 <- c('A', 'B', 'D', 'E')
obs_3 <- c('A', 'C', 'D', 'E')

f1 <- exp(forward(hmm, obs_1))
f2 <- exp(forward(hmm, obs_2))
f3 <- exp(forward(hmm, obs_3))

b1 <- exp(backward(hmm, obs_1))
b2 <- exp(backward(hmm, obs_2))
b3 <- exp(backward(hmm, obs_3))

p1 <- f1[-nrow(f1), -ncol(f1)] * b1[-nrow(b1), -ncol(b1)] / sum(f1[, ncol(f1)])
p2 <- f2[-nrow(f2), -ncol(f2)] * b2[-nrow(b2), -ncol(b2)] / sum(f2[, ncol(f2)])
p3 <- f3[-nrow(f3), -ncol(f3)] * b3[-nrow(b3), -ncol(b3)] / sum(f3[, ncol(f3)])

pc1 <- apply(p1, 1, sum)
pc2 <- apply(p2, 1, sum)
pc3 <- apply(p3, 1, sum)
sprintf("%1.2f", pc1)
sprintf("%1.2f", pc2)
sprintf("%1.2f", pc3)

pc <- pc1 + pc2 + pc3
sprintf("%1.2f", pc)