拓端tecdat|R語言中的隱馬爾可夫HMM模型實例

本文轉載自查看原文 2020-11-02 13:23 454 R語言/ HMM/ 隱馬爾可夫/ 廣義線性模型

原文鏈接：http://tecdat.cn/?p=17592

最近，我們使用隱馬爾可夫模型開發了一種解決方案，並被要求解釋這個方案。

HMM用於建模數據序列，無論是從連續概率分布還是從離散概率分布得出的。它們與狀態空間和高斯混合模型相關，因為它們旨在估計引起觀測的狀態。狀態是未知或“隱藏”的，並且HMM試圖估計狀態，類似於無監督聚類過程。

例子

在介紹HMM背后的基本理論之前，這里有一個示例，它將幫助您理解核心概念。有兩個骰子和一罐軟糖。B擲骰子，如果總數大於4，他會拿幾顆軟糖再擲一次。如果總數等於2，則他拿幾把軟糖，然后將骰子交給A。現在該輪到A擲骰子了。如果她的擲骰大於4，她會吃一些軟糖，但是她不喜歡黑色的其他顏色（兩極分化的看法），因此我們希望B會比A多。他們這樣做直到罐子空了。

現在假設A和B在不同的房間里，我們看不到誰在擲骰子。取而代之的是，我們只知道后來吃了多少軟糖。我們不知道顏色，僅是從罐子中取出的軟糖的最終數量。我們怎么知道誰擲骰子？HMM。

在此示例中，狀態是擲骰子的人，A或B。觀察結果是該回合中吃了多少軟糖。如果該值小於4，骰子的擲骰和通過骰子的條件就是轉移概率。由於我們組成了這個示例，我們可以准確地計算出轉移概率，即1/12。沒有條件說轉移概率必須相同，例如A擲骰子2時可以將骰子移交給他，例如，概率為1/36。

模擬

首先，我們將模擬該示例。B平均要吃12顆軟糖，而A則需要4顆。

# 設置
simulate <- function(N, dice.val = 6, jbns, switch.val = 4){
＃模擬變量
＃可以只使用一個骰子樣本
＃不同的機制，例如只丟1個骰子，或任何其他概率分布
b<- sample(1:dice.val, N, replace = T) + sample(1:dice.val, N, replace = T)
a <- sample(1:dice.val, N, replace = T) + sample(1:dice.val, N, replace = T)
bob.jbns <- rpois(N, jbns[1])
alice.jbns <- rpois(N, jbns[2])
# 狀態
draws <- data.frame(state = rep(NA, N), obs = rep(NA, N),
# 返回結果
return(cbind(roll = 1:N, draws))
# 模擬場景
draws <- simulate(N, jbns = c(12, 4), switch.val = 4)
# 觀察結果
ggplot(draws, aes(x = roll, y = obs)) + geom_line()

如您所見，僅檢查一系列計數來確定誰擲骰子是困難的。我們將擬合HMM。由於我們正在處理計數數據，因此觀察值是從泊松分布中得出的。

fit.hmm <- function(draws){
# HMM
mod <- fit(obs ~ 1, data = draws, nstates = 2, family = poisson()
# 通過估計后驗來預測狀態
est.states <- posterior(fit.mod)
head(est.states)
# 結果
hmm.post.df <- melt(est.states, measure.vars =
# 輸出表格
print(table(draws[,c("state", "est.state.labels")]))

## iteration 0 logLik: -346.2084
## iteration 5 logLik: -274.2033
## converged at iteration 7 with logLik: -274.2033
## est.state.labels
## state alice bob
## alice 49 2
## bob 3 46

模型迅速收斂。使用后驗概率，我們估計過程處於哪個狀態，即誰擁有骰子，A或B。要具體回答該問題，我們需要更多地了解該過程。在這種情況下，我們知道A只喜歡黑軟糖。否則，我們只能說該過程處於狀態1或2。下圖顯示了HMM很好地擬合了數據並估計了隱藏狀態。

# 繪圖輸出
g0 <- (ggplot(model.output$draws, aes(x = roll, y = obs)) + geom_line() +
theme(axis.ticks = element_blank(), axis.title.y = element_blank())) %>% ggplotGrob
g1 <- (ggplot(model.output$draws, aes(x = roll, y = state, fill = state, col = state)) +
g0$widths <- g1$widths
return(grid.arrange(g0, g1
plot.hmm.output(hmm1)

令人印象深刻的是，該模型擬合數據和濾除噪聲以估計狀態的良好程度。公平地說，可以通過忽略時間分量並使用EM算法來估計狀態。但是，由於我們知道數據形成一個序列，因為觀察下一次發生的概率取決於前一個即\（P（X_t | X_ {t-1}）\），其中\（X_t \ ）是軟糖的數量。

考慮到我們構造的問題，這可能是一個相對簡單的案例。如果轉移概率大得多怎么辦？

simulate(100, jbns = c(12, 4), switch.val = 7)

## iteration 0 logLik: -354.2707
## iteration 5 logLik: -282.4679
## iteration 10 logLik: -282.3879
## iteration 15 logLik: -282.3764
## iteration 20 logLik: -282.3748
## iteration 25 logLik: -282.3745
## converged at iteration 30 with logLik: -282.3745
## est.state.labels
## state alice bob
## alice 54 2
## bob 5 39

plot(hmm2)

這有很多噪音數據，但是HMM仍然做得很好。性能的提高部分歸因於我們對從罐中取出的軟糖數量的選擇。分布越明顯，模型就越容易拾取轉移。公平地講，我們可以計算中位數，並將所有低於中位數的值都歸為一個狀態，而將所有高於中位數的值歸為另一狀態，您可以從結果中看到它們做得很好。這是因為轉移概率非常高，並且預計我們會從每個狀態觀察到相似數量的觀察結果。當轉移概率不同時，我們會看到HMM表現更好。