公务员可以自己做网站吗/如何做营销策划方案

《统计学习方法》 李航著 第十章 隐马尔可夫模型

我是小白一个；本文代码转载地址文末有注释；有问题请多指教

先看书，看完书，代码就看懂了。程序只是将算法翻译成机器认识的罢了

import numpy as np

class HiddenMarkov:

def __init__(self):

self.alphas = None

self.forward_P = None

self.betas = None

self.backward_P = None

# 前向算法

# Q 是状态集合，里面包含了所有可能的状态

# V 是我们的观测的集合，里面包含了所有可能的观测结果

# A 状态转移概率分布

# B 观测概率分布

# O 观测序列，依次为观测值

# PI 初始概率分布。根据这个先生成初始状态。

def forward(self, Q, V, A, B, O, PI):

# 状态序列的大小

N = len(Q)

# 观测序列的大小

M = len(O)

# 初始化前向概率alpha值

alphas = np.zeros((N, M))

# 时刻数=观测序列数

T = M

# 遍历每一个时刻，计算前向概率alpha值

for t in range(T):

# 得到序列对应的索引

indexOfO = V.index(O[t])

# 遍历状态序列

for i in range(N):

# 初始化alpha初值

if t == 0:

# P176 公式(10.15)

alphas[i][t] = PI[t][i] * B[i][indexOfO]

print('alpha1(%d) = p%db%db(o1) = %f' %

(i + 1, i, i, alphas[i][t]))

else:

# P176 公式(10.16)

alphas[i][t] = np.dot([alpha[t - 1] for alpha in alphas],

[a[i] for a in A]) * B[i][indexOfO]

print('alpha%d(%d) = [sigma alpha%d(i)ai%d]b%d(o%d) = %f' %

(t + 1, i + 1, t - 1, i, i, t, alphas[i][t]))

# P176 公式(10.17)

self.forward_P = np.sum([alpha[M - 1] for alpha in alphas])

self.alphas = alphas

# 后向算法

# Q 是状态集合，里面包含了所有可能的状态

# V 是我们的观测的集合，里面包含了所有可能的观测结果

# A 状态转移概率分布

# B 观测概率分布

# O 观测序列，依次为观测值

# PI 初始概率分布。根据这个先生成初始状态。

def backward(self, Q, V, A, B, O, PI):

# 状态序列的大小

N = len(Q)

# 观测序列的大小

M = len(O)

# 初始化后向概率beta值，P178 公式(10.19)

betas = np.ones((N, M))

for i in range(N):

print('beta%d(%d) = 1' % (M, i + 1))

# 对观测序列逆向遍历

for t in range(M - 2, -1, -1):

# 得到序列对应的索引

indexOfO = V.index(O[t + 1])

# 遍历状态序列

for i in range(N):

# P178 公式(10.20)

betas[i][t] = np.dot(

np.multiply(A[i], [b[indexOfO] for b in B]),

[beta[t + 1] for beta in betas])

realT = t + 1

realI = i + 1

print('beta%d(%d) = sigma[a%djbj(o%d)beta%d(j)] = (' %

(realT, realI, realI, realT + 1, realT + 1),

end='')

for j in range(N):

print("%.2f * %.2f * %.2f + " %

(A[i][j], B[j][indexOfO], betas[j][t + 1]),

end='')

print("0) = %.3f" % betas[i][t])

# 取出第一个值

indexOfO = V.index(O[0])

self.betas = betas

# P178 公式(10.21)

P = np.dot(np.multiply(PI, [b[indexOfO] for b in B]),

[beta[0] for beta in betas])

self.backward_P = P

print("P(O|lambda) = ", end="")

for i in range(N):

print("%.1f * %.1f * %.5f + " %

(PI[0][i], B[i][indexOfO], betas[i][0]),

end="")

print("0 = %f" % P)

# 维特比算法：动态规划解隐马尔代夫模型预测问题

# Q 是状态集合，里面包含了所有可能的状态

# V 是我们的观测的集合，里面包含了所有可能的观测结果

# A 状态转移概率分布

# B 观测概率分布

# O 观测序列，依次为观测值

# PI 初始概率分布。根据这个先生成初始状态。

def viterbi(self, Q, V, A, B, O, PI):

# 状态序列的大小

N = len(Q)

# 观测序列的大小

M = len(O)

# 初始化daltas：存当前时刻当前状态的所有单个路径的概率最大值

deltas = np.zeros((N, M))

# 初始化psis：存当前时刻当前状态所有单个路径中概率最大路径的前一时刻结点

psis = np.zeros((N, M))

# 初始化最优路径矩阵，该矩阵维度与观测序列维度相同。这是我们最后的输出。

I = np.zeros((1, M))

# 遍历观测序列

for t in range(M):

# 递推从t=2开始

realT = t + 1

# 得到序列对应的索引

indexOfO = V.index(O[t])

for i in range(N):

realI = i + 1

if t == 0:

# P185 算法10.5 步骤(1)

deltas[i][t] = PI[0][i] * B[i][indexOfO]

psis[i][t] = 0

print('delta1(%d) = pi%d * b%d(o1) = %.2f * %.2f = %.2f' %

(realI, realI, realI, PI[0][i], B[i][indexOfO],

deltas[i][t]))

print('psis1(%d) = 0' % (realI))

else:

# # P185 算法10.5 步骤(2)

deltas[i][t] = np.max(

np.multiply([delta[t - 1] for delta in deltas],

[a[i] for a in A])) * B[i][indexOfO]

print(

'delta%d(%d) = max[delta%d(j)aj%d]b%d(o%d) = %.2f * %.2f = %.5f'

% (realT, realI, realT - 1, realI, realI, realT,

np.max(

np.multiply([delta[t - 1] for delta in deltas],

[a[i] for a in A])), B[i][indexOfO],

deltas[i][t]))

# 对于y=f(x)，argmax返回取得最大值y时的x

psis[i][t] = np.argmax(

np.multiply([delta[t - 1] for delta in deltas],

[a[i] for a in A]))

print('psis%d(%d) = argmax[delta%d(j)aj%d] = %d' %

(realT, realI, realT - 1, realI, psis[i][t]))

# 得到最优路径的终结点

I[0][M - 1] = np.argmax([delta[M - 1] for delta in deltas])

print('i%d = argmax[deltaT(i)] = %d' % (M, I[0][M - 1] + 1))

# 递归由后向前得到其他结点

for t in range(M - 2, -1, -1):

I[0][t] = psis[int(I[0][t + 1])][t + 1]

print('i%d = psis%d(i%d) = %d' %

(t + 1, t + 2, t + 2, I[0][t] + 1))

# 输出最优路径

print('最优路径是：', "->".join([str(int(i + 1)) for i in I[0]]))

# 习题10.1

Q = [1, 2, 3]

V = ['红', '白']

A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]

B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]

# O = ['红', '白', '红', '红', '白', '红', '白', '白']

O = ['红', '白', '红', '白'] # 习题10.1的例子

PI = [[0.2, 0.4, 0.4]]

HMM = HiddenMarkov()

HMM.forward(Q, V, A, B, O, PI)

print("P(O|λ)={}".format(HMM.forward_P))

# HMM.backward(Q, V, A, B, O, PI)

# print("P(O|λ)={}".format(HMM.backward_P))

# HMM.viterbi(Q, V, A, B, O, PI)

结果

前向算法

alpha1(1) = p0b0b(o1) = 0.100000

alpha1(2) = p1b1b(o1) = 0.160000

alpha1(3) = p2b2b(o1) = 0.280000

alpha2(1) = [sigma alpha0(i)ai0]b0(o1) = 0.077000

alpha2(2) = [sigma alpha0(i)ai1]b1(o1) = 0.110400

alpha2(3) = [sigma alpha0(i)ai2]b2(o1) = 0.060600

alpha3(1) = [sigma alpha1(i)ai0]b0(o2) = 0.041870

alpha3(2) = [sigma alpha1(i)ai1]b1(o2) = 0.035512

alpha3(3) = [sigma alpha1(i)ai2]b2(o2) = 0.052836

alpha4(1) = [sigma alpha2(i)ai0]b0(o3) = 0.021078

alpha4(2) = [sigma alpha2(i)ai1]b1(o3) = 0.025188

alpha4(3) = [sigma alpha2(i)ai2]b2(o3) = 0.013824

P(O|λ)=0.06009079999999999

后向算法

beta4(1) = 1

beta4(2) = 1

beta4(3) = 1

beta3(1) = sigma[a1jbj(o4)beta4(j)] = (0.50 * 0.50 * 1.00 + 0.20 * 0.60 * 1.00 + 0.30 * 0.30 * 1.00 + 0) = 0.460

beta3(2) = sigma[a2jbj(o4)beta4(j)] = (0.30 * 0.50 * 1.00 + 0.50 * 0.60 * 1.00 + 0.20 * 0.30 * 1.00 + 0) = 0.510

beta3(3) = sigma[a3jbj(o4)beta4(j)] = (0.20 * 0.50 * 1.00 + 0.30 * 0.60 * 1.00 + 0.50 * 0.30 * 1.00 + 0) = 0.430

beta2(1) = sigma[a1jbj(o3)beta3(j)] = (0.50 * 0.50 * 0.46 + 0.20 * 0.40 * 0.51 + 0.30 * 0.70 * 0.43 + 0) = 0.246

beta2(2) = sigma[a2jbj(o3)beta3(j)] = (0.30 * 0.50 * 0.46 + 0.50 * 0.40 * 0.51 + 0.20 * 0.70 * 0.43 + 0) = 0.231

beta2(3) = sigma[a3jbj(o3)beta3(j)] = (0.20 * 0.50 * 0.46 + 0.30 * 0.40 * 0.51 + 0.50 * 0.70 * 0.43 + 0) = 0.258

beta1(1) = sigma[a1jbj(o2)beta2(j)] = (0.50 * 0.50 * 0.25 + 0.20 * 0.60 * 0.23 + 0.30 * 0.30 * 0.26 + 0) = 0.112

beta1(2) = sigma[a2jbj(o2)beta2(j)] = (0.30 * 0.50 * 0.25 + 0.50 * 0.60 * 0.23 + 0.20 * 0.30 * 0.26 + 0) = 0.122

beta1(3) = sigma[a3jbj(o2)beta2(j)] = (0.20 * 0.50 * 0.25 + 0.30 * 0.60 * 0.23 + 0.50 * 0.30 * 0.26 + 0) = 0.105

P(O|lambda) = 0.2 * 0.5 * 0.11246 + 0.4 * 0.4 * 0.12174 + 0.4 * 0.7 * 0.10488 + 0 = 0.060091

P(O|λ)=[0.0600908]

维特比算法

delta1(1) = pi1 * b1(o1) = 0.20 * 0.50 = 0.10

psis1(1) = 0

delta1(2) = pi2 * b2(o1) = 0.40 * 0.40 = 0.16

psis1(2) = 0

delta1(3) = pi3 * b3(o1) = 0.40 * 0.70 = 0.28

psis1(3) = 0

delta2(1) = max[delta1(j)aj1]b1(o2) = 0.06 * 0.50 = 0.02800

psis2(1) = argmax[delta1(j)aj1] = 2

delta2(2) = max[delta1(j)aj2]b2(o2) = 0.08 * 0.60 = 0.05040

psis2(2) = argmax[delta1(j)aj2] = 2

delta2(3) = max[delta1(j)aj3]b3(o2) = 0.14 * 0.30 = 0.04200

psis2(3) = argmax[delta1(j)aj3] = 2

delta3(1) = max[delta2(j)aj1]b1(o3) = 0.02 * 0.50 = 0.00756

psis3(1) = argmax[delta2(j)aj1] = 1

delta3(2) = max[delta2(j)aj2]b2(o3) = 0.03 * 0.40 = 0.01008

psis3(2) = argmax[delta2(j)aj2] = 1

delta3(3) = max[delta2(j)aj3]b3(o3) = 0.02 * 0.70 = 0.01470

psis3(3) = argmax[delta2(j)aj3] = 2

delta4(1) = max[delta3(j)aj1]b1(o4) = 0.00 * 0.50 = 0.00189

psis4(1) = argmax[delta3(j)aj1] = 0

delta4(2) = max[delta3(j)aj2]b2(o4) = 0.01 * 0.60 = 0.00302

psis4(2) = argmax[delta3(j)aj2] = 1

delta4(3) = max[delta3(j)aj3]b3(o4) = 0.01 * 0.30 = 0.00220

psis4(3) = argmax[delta3(j)aj3] = 2

i4 = argmax[deltaT(i)] = 2

i3 = psis4(i4) = 2

i2 = psis3(i3) = 2

i1 = psis2(i2) = 3

最优路径是： 3->2->2->2