1. 程式人生 > 實用技巧 >R語言中的隱馬爾可夫HMM模型例項

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顆。

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

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

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

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

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

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

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

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

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

  1. ## iteration 0 logLik: -354.2707
  2. ## iteration 5 logLik: -282.4679
  3. ## iteration 10 logLik: -282.3879
  4. ## iteration 15 logLik: -282.3764
  5. ## iteration 20 logLik: -282.3748
  6. ## iteration 25 logLik: -282.3745
  7. ## converged at iteration 30 with logLik: -282.3745
  8. ## est.state.labels
  9. ## state alice bob
  10. ## alice 54 2
  11. ## bob 5 39

plot(hmm2)

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

如果觀察結果來自相同的分佈,即A和B吃了相同數量的軟糖怎麼辦?

  1. hmm3 <- fit.hmm(draws)
  2. plot(hmm3)

不太好,但這是可以預期的。如果從中得出觀察結果的分佈之間沒有差異,則可能也只有1個狀態。

實際如何估算狀態?

首先,狀態數量及其分佈方式本質上是未知的。利用對系統建模的知識,使用者可以選擇合理數量的狀態。在我們的示例中,我們知道有兩種狀態使事情變得容易。可能知道確切的狀態數,但這並不常見。再次通過系統知識來假設觀察結果通常是合理的,這通常是合理的。

從這裡開始,使用Baum-Welch演算法來估計引數,這是EM演算法的一種變體,它利用了觀測序列和Markov屬性。除了估計狀態的引數外,還需要估計轉移概率。Baum-Welch演算法首先對資料進行正向傳遞,然後進行反向傳遞。然後更新狀態轉移概率。然後重複此過程,直到收斂為止。

在現實世界

在現實世界中,HMM通常用於

  • 股票市場預測,無論市場處於牛市還是熊市
  • 估計NLP中的詞性
  • 生物測序
  • 序列分類

僅舉幾例。只要有觀察序列,就可以使用HMM,這對於離散情況也適用。


最受歡迎的見解

1.用R語言模擬混合制排隊隨機服務排隊系統

2.R語言中使用排隊論預測等待時間

3.R語言中實現馬爾可夫鏈蒙特卡羅MCMC模型

4.R語言中的馬爾科夫機制轉換(Markov regime switching)模型

5.matlab貝葉斯隱馬爾可夫hmm模型

6.用R語言模擬混合制排隊隨機服務排隊系統

7.Python基於粒子群優化的投資組合優化

8.R語言馬爾可夫轉換模型研究交通傷亡人數事故預測

9.用機器學習識別不斷變化的股市狀況——隱馬爾可夫模型的應用