#Hidden Markov Model dengan Pelatihan Menggunakan Algoritma Baum-Welch Untuk Maksimasi P(O|λ) Sebagai Dasar Pencarian Urutan Hidden State Terbaik Bersyarat Observasi Berganda Oleh Algoritma Viterbi
#data asli(misal)
set.seed(121450098)
test <- sample(1:2,5,replace=T)
test## [1] 1 1 2 2 2transisi <- matrix(c(0.5, 0.5, 0.4, 0.6), nrow = 2, byrow = TRUE)#ini matriks transisi antra state
p0 <- c(0.6, 0.4)#probabilitas inisiasi setiap state
matrix1 <- matrix(data=c(0.2,0.7,0.1,0.1,0.7,0.2), nrow = 2, byrow = TRUE)#matriks transisi hidden state ke variabel pengamatan ke-1
matrix2 <- matrix(data=c(0.1,0.8,0.05,0.05,0.25,0.5,0.05,0.2), nrow = 2, byrow = TRUE)#matriks transisi hidden state ke variabel pengamatan ke-2
matrix3 <- matrix(data=c(0.4,0.4,0.2,0.15,0.8,0.05), nrow = 2, byrow = TRUE)#matriks transisi hidden state ke variabel pengamatan ke-3
matrix4 <- matrix(data=c(0.2,0.7,0.1,0.1,0.7,0.2),nrow = 2,byrow = T)
ems <- list(matrix1, matrix2, matrix3,matrix4)#ini bua gabungin ketiga matriks tadi supaya lebih mudah
obs <- list(c(1,4,2,1), c(3,4,1,2), c(1,4,3,2),c(2,1,2,1),c(1,1,1,1))#ini observasinya
#Observasinya ada 3, setiap observasi melakukan pengamtan terhadap 3 variabel tersebut
#Dari setiap observasi dilihat mana hidden state yang paling mungkin terjadi pada waktu itu
#Penetapan suatu state pada suatu observasi pada waktu ke-t tentunya juga memengaruhi penetapan state di waktu ke t+1. Oleh karena itu, diperlukan peninjauan terhadap urutan hidden state yang tepat dari waktu ke-1 sampai ke T(akhir). Cara melakukannya adalah dengan mencari urutan hidden state yang mampu memkasimalkan P(HiddenState_1:T | Observasi_1:T). Menggunakan cara brute force adalah salah satu cara sederhananya. Namun, beban kemangkusan menjadi taruhannya. Brute force melakuakn percobaan sebanyak observasi yakni sebanyak n_state^T. Membandingkan satu-persatu P(HiddenState_1:T | Observasi_1:T).
ems## [[1]]
## [,1] [,2] [,3]
## [1,] 0.2 0.7 0.1
## [2,] 0.1 0.7 0.2
##
## [[2]]
## [,1] [,2] [,3] [,4]
## [1,] 0.10 0.8 0.05 0.05
## [2,] 0.25 0.5 0.05 0.20
##
## [[3]]
## [,1] [,2] [,3]
## [1,] 0.40 0.4 0.20
## [2,] 0.15 0.8 0.05
##
## [[4]]
## [,1] [,2] [,3]
## [1,] 0.2 0.7 0.1
## [2,] 0.1 0.7 0.2viterbimultobs <- function(transisi,emisi,p0,seq_observasi){
pengamatan <- length(seq_observasi)
state <- length(p0)
deltas <- matrix(0,nrow=pengamatan,ncol=state)
psi <- matrix(0,nrow=pengamatan,ncol=state)
pmult <- function(t,obs,j,ems){
kali <- 1
for (i in 1:length(obs[[t]])) {
kali <- kali * ems[[i]][j,obs[[t]][i]]
}
return(kali)
}
for(i in 1:state){
deltas[1,i]<-pmult(1,obs,i,ems)
}
for(t in 2:pengamatan){
for(j in 1:state){
maxvector1 <- c()
for (i in 1:state) {
cekmax1 <- transisi[i,j] * deltas[t-1,i]
maxvector1<-append(maxvector1,cekmax1)
}
deltas[t,j] <- max(maxvector1)*pmult(t,seq_observasi,j,emisi)
}
}
for (t in 2:pengamatan) {
for(j in 1:state){
maxvector2<-c()
for(i in 1:state){
cekmax2 <- transisi[i,j] * deltas[t-1,i]
maxvector2 <- append(maxvector2,cekmax2)
}
psi[t,j] <- which.max(maxvector2)
}
}
qstar <- numeric(pengamatan)
qstar[pengamatan] <- which.max(deltas[pengamatan, ])
for (i in (pengamatan - 1):1) {
qstar[i] <- psi[i + 1, qstar[i + 1]]
}
return(qstar)
}
viterbimultobs(transisi, ems, p0, obs)## [1] 2 2 1 2 1##Forward dan Backward
forwardmultobs<- function(transisi,emisi,p0,seq_observasi) {
pengamatan <- length(seq_observasi)
state <- length(p0)
alpha<-matrix(0,pengamatan,state)
pmult <- function(t,obs,j,ems){
kali <- 1
for (i in 1:length(obs[[t]])) {
kali <- kali * ems[[i]][j,obs[[t]][i]]
}
return(kali)
}
for(i in 1:state){
alpha[1,i]<-pmult(1,obs,i,ems) * p0[i]
}
for (t in 2:pengamatan) {
for (j in 1:state) {
alpha[t,j]<-0
for (i in 1:state){
alpha[t,j] <- alpha[t,j] + transisi[i,j] * alpha[t-1,i] * pmult(t,seq_observasi,j,emisi)
}
}
}
pstar <- sum(alpha[pengamatan,])
return(list(pstar = pstar, alpha=alpha))
}
forwardmultobs(transisi,ems,p0,obs)## $pstar
## [1] 2.996819e-14
##
## $alpha
## [,1] [,2]
## [1,] 4.800000e-04 6.400000e-04
## [2,] 6.944000e-07 2.620800e-06
## [3,] 1.953728e-09 1.343776e-09
## [4,] 8.480497e-12 2.496381e-11
## [5,] 2.276124e-14 7.206951e-15backwardmultobs<- function(transisi,emisi,p0,seq_observasi) {
state <- length(p0)
beta <- matrix(0,length(seq_observasi),state)
beta[length(seq_observasi), ] <- 1
pmult <- function(t,obs,j,ems){
kali <- 1
for (i in 1:length(obs[[t]])) {
kali <- kali * ems[[i]][j,obs[[t]][i]]
}
return(kali)
}
for (t in (length(seq_observasi)- 1):1) {
for (i in 1:state) {
beta[t,i] <- 0
for (j in 1:state) {
beta[t, i] <- beta[t,i] + transisi[i,j] * pmult(t+1,seq_observasi,j,emisi) * beta[t + 1,j]
}
}
}
pstar <- 0
for (i in 1:state) {
pstar <- pstar + p0[i] * pmult(1, seq_observasi, i, emisi) * beta[1, i]
}
return(list(pstar = pstar, beta=beta))
}
backwardmultobs(transisi,ems,p0,obs)## $pstar
## [1] 2.996819e-14
##
## $beta
## [,1] [,2]
## [1,] 2.537583e-11 2.779343e-11
## [2,] 9.491300e-09 8.919960e-09
## [3,] 8.820000e-06 9.478000e-06
## [4,] 9.875000e-04 8.650000e-04
## [5,] 1.000000e+00 1.000000e+00#Fungsinya untuk melatih sampai sesuai yang diinginkan
em_pstar_with_tolerance <- function(transisi, ems, p0, seq_observasi) {
pstar <- forwardmultobs(transisi, ems, p0, seq_observasi)$pstar
pmult <- function(t,obs,j,ems){
kali <- 1
for (i in 1:length(obs[[t]])) {
kali <- kali * ems[[i]][j,obs[[t]][i]]
}
return(kali)
}
state <- length(p0)
pengamatan <- length(seq_observasi)
# simpan dulu hasil alpha dan beta selayaknya sudah menjalankan fase forward dan backward
alpha <- forwardmultobs(transisi, ems, p0, seq_observasi)$alpha
beta <- backwardmultobs(transisi, ems, p0, seq_observasi)$beta
# fase perhitungan ksi dan gamma
# mulai fase perhitungan ksi
ksi_denom <- matrix(0,nrow=pengamatan-1,ncol=1)
for (t in 1:(pengamatan - 1)) {
for (i in 1:state) {
for (j in 1:state) {
ksi_denom[t] <- ksi_denom[t] + (alpha[t, i] * transisi[i, j] * pmult(t+1,obs,j,ems) * beta[t + 1, j])
}
}
}
ksi <- array(0, dim = c(pengamatan - 1, state, state))
for (t in 1:(pengamatan - 1)) {
for (i in 1:state) {
for (j in 1:state) {
ksi[t, i, j] <- (alpha[t, i] * transisi[i, j] * pmult(t+1,obs,j,ems) * beta[t + 1, j]) / ksi_denom[t]
}
}
}
# mulai fase perhitungan gamma
gamma <- matrix(0, nrow = pengamatan, ncol = state)
for (t in 1:pengamatan) {
denom_gamma <- 0
for (i in 1:state) {
denom_gamma <- denom_gamma + (alpha[t,i ] * beta[t, i])
}
for (i in 1:state) {
gamma[t, i] <- (alpha[t, i] * beta[t, i]) / denom_gamma
}
}
# Mulai fase pembaharuan transisi, emisi, dan p0
# fase pembaharuan p0
p0_topi <- gamma[1, ]
# fase pembaharuan transisi (harus dinormalisasi tentunya)
transisi_topi <- matrix(0, nrow = state, ncol = state)
for (i in 1:state) {
denom_transisi_topi <- sum(gamma[1:pengamatan-1, i])
for (j in 1:state) {
for (t in 1:(pengamatan - 1)) {
transisi_topi[i, j] <- transisi_topi[i, j] + (ksi[t, i, j] / denom_transisi_topi)
}
}
}
for (i in 1:state) { # normalisasi dulu
normalize_markov <- sum(transisi_topi[i, ])
for (j in 1:state) {
transisi_topi[i, j] <- transisi_topi[i, j] / normalize_markov
}
}
# fase pembaharuan emisi
delta <- function(t,k,v,seq_observasi) { # fungsi delta, kalau observasi ke-t != k maka akan return 0
if (seq_observasi[[t]][k] == v) {
return(1)
} else {
return(0)
}
}
normalize <- function(matriks) {
for (i in 1:nrow(matriks)) {
denom_normalize <- sum(matriks[i, ])
if (denom_normalize > 0) {
matriks[i, ] <- matriks[i, ] / denom_normalize
}
}
return(matriks)
}
dimensi <- length(seq_observasi[[1]])
list_emisi <- vector("list", dimensi)
for (k in 1:dimensi) {
emisi_topi <- matrix(0, nrow = state, ncol = ncol(ems[[k]]))
for (j in 1:state) {
denom_emisi_topi <- sum(gamma[, j])
for (v in 1:ncol(ems[[k]])) {
for (t in 1:pengamatan) {
emisi_topi[j, v] <- emisi_topi[j, v] + (gamma[t, j] * delta(t,k,v,seq_observasi)) / denom_emisi_topi
}
}
}
list_emisi[[k]] <- normalize(emisi_topi)
}
#untuk hitung pstar baru
pstar_topi <- forwardmultobs(transisi_topi,list_emisi,p0_topi,seq_observasi)$pstar
return(list(p0_new = p0_topi,
emisi_new = list_emisi,
transisi_new = transisi_topi,
pstar_topi = pstar_topi))
}
result<-em_pstar_with_tolerance(transisi,ems,p0,obs)#memunculkan parameter yang sudah dilatih dan pstar yang meningkat
afterEM <- viterbimultobs(result$transisi_new,result$emisi_new,result$p0_new,obs)
cat("Setelah EM urutannya :",afterEM)#kasus ini sama kebetulan baik sebelum dengan sesudah ## Setelah EM urutannya : 2 2 1 2 1print("Hasil sesudah EM")## [1] "Hasil sesudah EM"result#setelah pelatihan(di dalamnya mencakup pstar dan juga parameter yang dilatih)## $p0_new
## [1] 0.4064442 0.5935558
##
## $emisi_new
## $emisi_new[[1]]
## [,1] [,2] [,3]
## [1,] 0.7770997 0.1247341 0.09816622
## [2,] 0.4562280 0.2611019 0.28267006
##
## $emisi_new[[2]]
## [,1] [,2] [,3] [,4]
## [1,] 0.4637519 0 0 0.5362481
## [2,] 0.3482453 0 0 0.6517547
##
## $emisi_new[[3]]
## [,1] [,2] [,3]
## [1,] 0.4371841 0.3061553 0.2566607
## [2,] 0.3698135 0.4761844 0.1540021
##
## $emisi_new[[4]]
## [,1] [,2] [,3]
## [1,] 0.6451731 0.3548269 0
## [2,] 0.5633279 0.4366721 0
##
##
## $transisi_new
## [,1] [,2]
## [1,] 0.4430792 0.5569208
## [2,] 0.4675211 0.5324789
##
## $pstar_topi
## [1] 8.559356e-10print("Hasil sebelum EM")## [1] "Hasil sebelum EM"forwardmultobs(transisi,ems,p0,obs)$pstar## [1] 2.996819e-14list(transisi,ems,p0)#sebelum pelatihan## [[1]]
## [,1] [,2]
## [1,] 0.5 0.5
## [2,] 0.4 0.6
##
## [[2]]
## [[2]][[1]]
## [,1] [,2] [,3]
## [1,] 0.2 0.7 0.1
## [2,] 0.1 0.7 0.2
##
## [[2]][[2]]
## [,1] [,2] [,3] [,4]
## [1,] 0.10 0.8 0.05 0.05
## [2,] 0.25 0.5 0.05 0.20
##
## [[2]][[3]]
## [,1] [,2] [,3]
## [1,] 0.40 0.4 0.20
## [2,] 0.15 0.8 0.05
##
## [[2]][[4]]
## [,1] [,2] [,3]
## [1,] 0.2 0.7 0.1
## [2,] 0.1 0.7 0.2
##
##
## [[3]]
## [1] 0.6 0.4