统计学习方法第十章课后习题10.1,10.2,10.3

代码经书上正文中例子测试，通过

10.1本题要求用后向计算，这里同时提供正向算法代码

• import numpy as np
  
  def back_prob():
      beta = np.ones(N)
      for t in range(T - 2, -1, -1):
          beta_old = beta.copy()
          for i in range(N):
              beta[i] = 0
              for j in range(N):
                  beta[i] += A[i, j] * B[j, O[t + 1] ] * beta_old[j]
      p = 0
      for i in range(N):
          p += PI[i] * B[i, O[0]] * beta[i]
      return p
  
  def forward_prob():
      alpha = PI*B[:,O[0]]
      for t in range(T-1):
          alpha_old = alpha.copy()
          for i in range(N):
              alpha[i] = 0
              for j in range(N):
                  alpha[i] += alpha_old[j]*A[j,i]
              alpha[i] *= B[ i,O[t+1] ]
      p = sum(alpha)
      return p
  
  A = np.array([
      [0.5, 0.2, 0.3],
      [0.3, 0.5, 0.2],
      [0.2, 0.3, 0.5]
  ])
  
  B = np.array([
      [0.5, 0.5],
      [0.4, 0.6],
      [0.7, 0.3]
  ])
  
  PI = np.array([0.2, 0.4, 0.4])
  O = [0, 1, 0, 1]  # 0表示红，1表示白
  T = 4
  N = 3
  M = 2
  
  p = back_prob()
  print("后向 P(O|lambda) = ", p)
  p = forward_prob()
  print("前向 P(O|lambda) = ", p)
  

10.2  这里的概率，其实就是γ_t4(q3) 

按照书上公式（10.24），前向后向过程记录alpha,beta

import numpy as np

def back():
    BETA = []
    beta = np.ones(N)
    BETA.append(beta)
    for t in range(T - 2, -1, -1):
        beta_old = beta.copy()
        for i in range(N):
            beta[i] = 0
            for j in range(N):
                beta[i] += A[i, j] * B[j, O[t + 1] ] * beta_old[j]
        BETA.append(beta)

    return np.array(BETA).reshape(T, N)


def forward():
    ALPHA = []
    alpha = PI*B[:,O[0]]
    ALPHA.append(alpha)
    for t in range(T-1):
        alpha_old = alpha.copy()
        for i in range(N):
            alpha[i] = 0
            for j in range(N):
                alpha[i] += alpha_old[j]*A[j,i]
            alpha[i] *= B[ i,O[t+1] ]
        ALPHA.append(alpha)

    return np.array(ALPHA).reshape(T,N)


A = np.array([
    [0.5,0.1,0.4],
    [0.3,0.5,0.2],
    [0.2,0.2,0.6]
    ])
B = np.array([
    [0.5,0.5],
    [0.4,0.6],
    [0.7,0.3]
    ])
PI = np.array([0.2,0.3,0.5])
T = 8
N = 3
M = 2
O = np.array([0,1,0,0,1,0,1,1])
BETA = back()
ALPHA = forward()
p = ALPHA[3,2]*BETA[3,2] / sum(ALPHA[3,:]*BETA[3,:])
print("P(i4=q3|O,lambda) = ",p)

10.3 

最优路径：

import numpy as np

def viterbi():
    delta = PI * B[:, O[0]]
    FAI = []
    fai = np.zeros(N)
    FAI.append(fai)
    for t in range(1, T):
        delta_old = delta.copy()
        for i in range(N):
            delta[i] = np.max(delta_old * A[:, i] * B[i, O[t]])
            fai[i] = np.argmax(delta_old * A[:, i])
        FAI.append(fai.copy())  #注意，应为copy
    FAI = np.array(FAI).reshape(T, N)
    path = []
    I = np.argmax(delta)
    path.append(I)
    for t in range(T - 2, -1, -1):
        I = int( FAI[t + 1, I] )
        path.append(I)
    return (np.array(path)+1)[::-1] # 使状态从1开始计数 状态正向顺序


A = np.array([
    [0.5, 0.2, 0.3],
    [0.3, 0.5, 0.2],
    [0.2, 0.3, 0.5]
])

B = np.array([
    [0.5, 0.5],
    [0.4, 0.6],
    [0.7, 0.3]
])

PI = np.array([0.2, 0.4, 0.4])
O = np.array([0, 1, 0, 1])  # 0表示红，1表示白
T = 4
# O = np.array([0, 1, 0])  # 0表示红，1表示白
# T = 3
N = 3
M = 2

path = viterbi()
print(path)