本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程
你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】
01 隐马尔可夫模型
HMM(Hidden Markov Model), 也称隐性马尔可夫模型，是一个概率模型，用来描述一个系统隐性状态的转移和隐性状态的表现概率。
既是马尔可夫模型，就一定存在马尔可夫链，该马尔可夫链服从马尔可夫性质：即无记忆性。也就是说，这一时刻的状态，受且只受前一时刻的影响，而不受更往前时刻的状态的影响。
系统的隐性状态指的就是一些外界不便观察(或观察不到)的状态,
比如在当前的例子里面, 系统的状态指的是大叔使用骰子的状态，即{正常骰子, 作弊骰子1, 作弊骰子2,...}隐性状态的表现也就是,
可以观察到的，由隐性状态产生的外在表现特点。这里就是说, 骰子掷出的点数. {1,2,3,4,5,6}HMM模型将会描述，系统隐性状态的转移概率。也就是大叔切换骰子的概率,下图是一个例子，这时候大叔切换骰子的可能性被描述得淋漓尽致。
02 隐马尔可夫模型的解释
最经典的例子，掷骰子。假设我手里有三个不同的骰子。第一个骰子是我们平常见的骰子（称这个骰子为D6），6个面，每个面（1，2，3，4，5，6）出现的概率是1/6。第二个骰子是个四面体（称这个骰子为D4），每个面（1，2，3，4）出现的概率是1/4。第三个骰子有八个面（称这个骰子为D8），每个面（1，2，3，4，5，6，7，8）出现的概率是1/8。
假设我们开始掷骰子，我们先从三个骰子里挑一个，挑到每一个骰子的概率都是1/3。然后我们掷骰子，得到一个数字，1，2，3，4，5，6，7，8中的一个。不停的重复上述过程，我们会得到一串数字，每个数字都是1，2，3，4，5，6，7，8中的一个。例如我们可能得到这么一串数字（掷骰子10次）：1 6 3 5 2 7 3 5 2 4这串数字叫做可见状态链。但是在隐马尔可夫模型中，我们不仅仅有这么一串可见状态链，还有一串隐含状态链。在这个例子里，这串隐含状态链就是你用的骰子的序列。比如，隐含状态链有可能是：D6 D8 D8 D6 D4 D8 D6 D6 D4 D8一般来说，HMM中说到的马尔可夫链其实是指隐含状态链，因为隐含状态（骰子）之间存在转换概率（transition probability）。在我们这个例子里，D6的下一个状态是D4，D6，D8的概率都是1/3。D4，D8的下一个状态是D4，D6，D8的转换概率也都一样是1/3。这样设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，我们可以这样定义，D6后面不能接D4，D6后面是D6的概率是0.9，是D8的概率是0.1。
这样就是一个新的HMM。同样的，尽管可见状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率（emission probability）。就我们的例子来说，六面骰（D6）产生1的输出概率是1/6。产生2，3，4，5，6的概率也都是1/6。我们同样可以对输出概率进行其他定义。比如，我有一个被赌场动过手脚的六面骰子，掷出来是1的概率更大，是1/2，掷出来是2，3，4，5，6的概率是1/10。
03 隐马尔科夫模型R语言实现
## 隐马尔可夫模型的简单应用
## 状态序列:隐藏马尔科夫链随机生成的状态序列，称为i状态序列
## 观测序列：每个状态生成一个观测，由此产生的观测是随机的序列，称为观测序列
## 隐马尔科夫模型由，初始概率分布、状态转移概率分布以及观测概率分布确定
## 问题：我们有三个盒子，每个盒子都放有红球和白球
## 问题来自《统计学习方法》
## HMM - Hidden Markov Models
library(HMM)
## 所有可能状态集合
state <- c("1","2","3")
## 所有可能观测集合
Symbols <- c("红","白")
## 初始概率分布：第一次抽取到相应盒子的概率
startprob <- c(0.2,0.4,0.4)
## 状态转移概率矩阵
transpro <- matrix(c(0.5, 0.3, 0.2, 0.2, 0.5, 0.3,0.3,0.2,0.5),nrow = 3)
transpro
## [,1] [,2] [,3]
## [1,] 0.5 0.2 0.3
## [2,] 0.3 0.5 0.2
## [3,] 0.2 0.3 0.5
## 观测概率分布：每个盒子里不同颜色球的概率
emisspr <- matrix(c(0.5,0.4,0.7,0.5,0.6,0.3),nrow = 3)
emisspr
## [,1] [,2]
## [1,] 0.5 0.5
## [2,] 0.4 0.6
## [3,] 0.7 0.3
## 初始化隐马尔科夫模型
hmm <- initHMM(States = state,Symbols = Symbols,startProbs = startprob,transProbs = transpro,
emissionProbs = emisspr)
## 的到的一个观测
ober <- c("红","白","红","红")
## 1:概率计算问题
## 使用向前算法的到相应观测的概率
## 包含向前概率的矩阵。第一个维度指的是状态和第二维度的时间。
exp(forward(hmm, ober))
## index
## states 1 2 3 4
## 1 0.10 0.0770 0.041870 0.02107790
## 2 0.16 0.1104 0.035512 0.01679232
## 3 0.28 0.0606 0.052836 0.03225698
## 状态求和，得到每次个时间为该序列的概率
colSums(exp(forward(hmm, ober)))
## 1 2 3 4
## 0.5400000 0.2480000 0.1302180 0.0701272
## 使用向后算法的到相应观测的概率
## 包含向后概率的矩阵。第一个维度指的是状态和第二维度的时间。
exp(backward(hmm, ober))
## index
## states 1 2 3 4
## 1 0.132638 0.2939 0.54 1
## 2 0.140463 0.2588 0.49 1
## 3 0.122819 0.3123 0.57 1
colSums(exp(backward(hmm, ober)))
## 1 2 3 4
## 0.39592 0.86500 1.60000 3.00000
## 预测问题：给定观测序列，求最有可能的对应的状态序列
## 维特比算法求解
## 输出：一个字符串向量，包含最可能的状态路径。
viterbi(hmm,ober)
## [1] "3" "3" "3" "3"
## 此函数计算在一个给定的序列的观测和给定的Hidden Markov模型在时间k在状态x的后验概率。
posterior(hmm,ober)
## index
## states 1 2 3 4
## 1 0.1891392 0.3227036 0.3224113 0.3005667
## 2 0.3204759 0.4074242 0.2481331 0.2394552
## 3 0.4903849 0.2698722 0.4294556 0.4599782
## 从结果可以得出，最有可能的状态序列为：3，2，3，3
colSums(posterior(hmm,ober))
## 1 2 3 4
## 1 1 1 1
## 给定一个隐马尔科夫模型的模拟
simHMM(hmm,4)
## $states
## [1] "3" "2" "2" "1"
##
## $observation
## [1] "红" "白" "白" "红"
## 学习问题：已知观测序列，学习模型的参数
hmm <- initHMM(States = state,Symbols = Symbols,startProbs = startprob,
transProbs = transpro,emissionProbs = emisspr)
hmm
## $States
## [1] "1" "2" "3"
##
## $Symbols
## [1] "红" "白"
##
## $startProbs
## 1 2 3
## 0.2 0.4 0.4
##
## $transProbs
## to
## from 1 2 3
## 1 0.5 0.2 0.3
## 2 0.3 0.5 0.2
## 3 0.2 0.3 0.5
##
## $emissionProbs
## symbols
## states 红 白
## 1 0.5 0.5
## 2 0.4 0.6
## 3 0.7 0.3
ober <- c(sample(c(rep("红",100),rep("白",200))),
sample(c(rep("白",300),rep("红",100))))
bw <- baumWelch(hmm,observation = ober)
bw$hmm
## $States
## [1] "1" "2" "3"
##
## $Symbols
## [1] "红" "白"
##
## $startProbs
## 1 2 3
## 0.2 0.4 0.4
##
## $transProbs
## to
## from 1 2 3
## 1 0.5162805 0.2293114 0.2544081
## 2 0.2987254 0.5283639 0.1729107
## 3 0.2226074 0.4078684 0.3695243
##
## $emissionProbs
## symbols
## states 红 白
## 1 0.2942211 0.7057789
## 2 0.2399061 0.7600939
## 3 0.3448797 0.6551203
## 连续迭代的差异
bw$difference
## [1] 0.711033336 0.006630626 0.006323676 0.006074176 0.005842350
## [6] 0.005626430 0.005424877 0.005236343 0.005059639 0.004893716
## [11] 0.004737642 0.004590589 0.004451817 0.004320666 0.004196539
## [16] 0.004078903 0.003967275 0.003861218 0.003760336 0.003664268
## [21] 0.003572688 0.003485297 0.003401820 0.003322007 0.003245628
## [26] 0.003172473 0.003102346 0.003035068 0.002970474 0.002908409
## [31] 0.002848733 0.002791312 0.002736026 0.002682760 0.002631408
## [36] 0.002581872 0.002534060 0.002487884 0.002443265 0.002400128
## [41] 0.002358401 0.002318018 0.002278918 0.002241040 0.002204332
## [46] 0.002168739 0.002134215 0.002100712 0.002068187 0.002036598
## [51] 0.002005908 0.001976079 0.001947077 0.001918868 0.001891421
## [56] 0.001864706 0.001838696 0.001813363 0.001788683 0.001764630
## [61] 0.001741182 0.001718317 0.001696014 0.001674253 0.001653014
## [66] 0.001632281 0.001612035 0.001592261 0.001572941 0.001554062
## [71] 0.001535608 0.001517565 0.001499921 0.001482663 0.001465778
## [76] 0.001449255 0.001433083 0.001417251 0.001401748 0.001386564
## [81] 0.001371691 0.001357119 0.001342838 0.001328841 0.001315120
## [86] 0.001301666 0.001288472 0.001275531 0.001262836 0.001250379
## [91] 0.001238155 0.001226157 0.001214379 0.001202815 0.001191460
## [96] 0.001180307 0.001169353 0.001158591 0.001148016 0.001137625
参考资料
https://www.zhihu.com/question/20962240
关注R小盐，关注科研私家菜(VX_GZH: SciPrivate)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型
