隨機條件場—概率計算問題、IIS/GD學習算法、維特比預測算法

這一章的算法不是很好寫，整整研究了好幾天，代碼還是有點小問題，僅供參考。
用的是書上定義的特征函數。

import numpy as np class CRF:def __init__(self,y=None,x=None,y_num=None,x_num=None,N=None):self.y = yself.x = xself.y_num = y_numself.x_num = x_numself.N = Nself.get_feature()self.build_Marix(self.x[0])def get_feature(self):self.ti = [lambda y_1, y, x, i: 1 if i == 2 and y_1 == 1 and y == 2 else 0,lambda y_1, y, x, i: 1 if i == 3 and y_1 == 1 and y == 2 else 0,lambda y_1, y, x, i: 1 if i == 2 and y_1 == 1 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 3 and y_1 == 2 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 2 and y_1 == 2 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 3 and y_1 == 2 and y == 2 else 0,]self.w_ti = [1,1,0.6,1,1,0.2]self.si = [lambda y_1, y, x, i: 1 if i == 1 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 1 and y == 2 else 0,lambda y_1, y, x, i: 1 if i == 2 and y == 2 else 0,lambda y_1, y, x, i: 1 if i == 2 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 3 and y == 1 else 0,lambda y_1, y, x, i: 1 if i == 3 and y == 2 else 0,]self.w_si = [1,0.5,0.5,0.8,0.8,0.5]self.fk = self.ti+self.siself.wk = self.w_ti+self.w_sidef build_Marix(self,x):self.Marix = np.zeros((self.N+1,self.y_num,self.y_num))for i in range(self.N+1):for n in range(self.y_num):if i == self.N:self.Marix[i][:,0] = 1breakfor m in range(self.y_num):for k in range(len(self.fk)):if i == 0:if n==0:self.Marix[i][0][m] += self.wk[k] * (self.fk[k](0,m+1,x[i],i+1))else:self.Marix[i][n][m] += self.wk[k] * (self.fk[k](n+1,m+1,x[i],i+1))self.Marix = np.exp(self.Marix)def get_aiT(self):self.aiT = np.zeros((self.N+1,self.y_num))self.aiT[0,:] = 1for i in range(1,self.N+1):self.aiT[i] = self.aiT[i-1].dot(self.Marix[i])def get_biT(self):self.biT = np.zeros((self.N + 1, self.y_num))self.biT[self.N,:] = 1for i in range(self.N-1,-1,-1):self.biT[i] = self.Marix[i+1].dot(self.biT[i+1])#--------概率計算問題---------------def predict_Py_i(self,i,yi):self.get_aiT()self.get_biT()return self.aiT[i][yi-1]*self.biT[i][yi-1]/np.sum(self.aiT[self.N])def compute_p(self,y):result = 1for i in range(self.N):if i == 0:Z = self.Marix[i]result *= self.Marix[i][0][y[i]-1]else:Z = Z.dot(self.Marix[i])result *= self.Marix[i][y[i-1]-1][y[i]-1]return result/np.sum(Z)#----------維特比預測算法-------------------------def Viterbi(self,x):p_marix = np.zeros((len(x),self.y_num))rout = [[0] for i in range(self.y_num)]for j in range(self.y_num):p_marix[0][j] = np.sum(np.array([f(0,j+1,x[0],1) for f in self.fk])*self.wk)for i in range(1,len(x)):for l in range(self.y_num):max = 0rout_tem = 0for j in range(self.y_num):wf = p_marix[i - 1][j] + np.sum(np.array([f(j+1, l+1, x[i], i + 1) for f in self.fk]) * self.wk)if wf > max:max = wfrout_tem = jp_marix[i][l] = maxrout[l].append(rout_tem)max_result = np.max(p_marix[len(x)-1])max_index = list(p_marix[len(x)-1]).index(np.max(p_marix[len(x)-1]))rout[max_index].append(max_index)return max_result,[i+1 for i in rout[max_index][1:]]# -----------學習算法--------------------------------def compute_pxy_f(self):#pfklist_pxyti = np.zeros(len(self.ti))list_pxysi = np.zeros(len(self.si))#p_fklist_countti = np.zeros(len(self.ti))list_countsi = np.zeros(len(self.si))for x,y in zip(self.x,self.y):self.build_Marix(x)self.get_aiT()self.get_biT()for i in range(len(x)):list_temp = np.zeros(len(self.ti))for k in range(len(self.ti)):if i == 0:right = 1*self.Marix[0][0][y[i]-1]*self.biT[0][y[i]-1]left = self.ti[k](i,y[i],x[i],i+1)list_temp[k] += left*right/np.sum(self.aiT[self.N])else:right = self.aiT[i][y[i-1]-1] * self.Marix[i][y[i-1]-1][y[i]-1]*self.biT[i][y[i]-1]left = self.ti[k](y[i-1],y[i],x[i],i+1)list_temp[k] += left*right/np.sum(self.aiT[self.N])if left == 1:list_countti[k] += 1list_pxyti += list_tempfor i in range(len(x)):list_temp = np.zeros(len(self.si))for k in range(len(self.si)):right = self.aiT[i][y[i]-1]*self.biT[i][y[i]-1]left = self.si[k](y[i-1], y[i], x[i], i+1)list_temp[k] += left * right / np.sum(self.aiT[self.N])if left == 1:list_countsi[k] += 1list_pxysi += list_tempreturn list_countti/len(self.x),list_countsi/len(self.x),list_pxyti/len(self.x),list_pxysi/len(self.x)def compute_fw(self):left = 0right = 0for x in self.x:self.build_Marix(x)self.get_aiT()left += np.log(np.sum(self.aiT[self.N]))for x, y in zip(self.x, self.y):for i in range(len(x)):for k in range(len(self.fk)):if i == 0:right += self.wk[k] * self.fk[k](0,y[i],x[i],i+1)else:right += self.wk[k] * self.fk[k](y[i-1], y[i], x[i], i+1)return (left - right)/len(self.x)def fit(self,max_iter=3,how='IIS',lr=0.001):self.w_ti = [0]*len(self.ti)self.w_si = [0]*len(self.si)self.wk = self.w_ti + self.w_siself.x = np.array(self.x)self.y = np.array(self.y)if how == 'IIS':S = 20for i in range(max_iter):ep_tk,ep_sk,eptk,epsk = self.compute_pxy_f()if np.linalg.norm(1/S*np.log(ep_tk/eptk), ord=2)+np.linalg.norm(1/S*np.log(ep_sk/epsk),ord=2) < 0.1:print('when iter is '+str(i)+' shoulian')breakself.w_ti += 1/S*np.log(ep_tk/eptk)self.w_si += 1/S*np.log(ep_sk/epsk)self.wk = list(self.w_ti) + list(self.w_si)#按著最大熵模型的公式寫的，用EP - E_P為倒數，不知道對不對elif how == 'GD':for i in range(max_iter):ep_tk,ep_sk,eptk,epsk = self.compute_pxy_f()gt = eptk - ep_tkgs = epsk - ep_skfanshut = np.linalg.norm(gt,ord=2)fanshus = np.linalg.norm(gs,ord=2)if fanshut + fanshus < 0.75:print('when iter is '+str(i)+' shoulian')breaktemp_w_ti = self.w_titemp_w_si = self.w_sifw_list = []#線性搜索這里有問題，fw會一直縮小，因該是fw的計算出現了錯誤for k in range(20):self.w_ti = temp_w_ti - lr * gt * kself.w_si = temp_w_si - lr * gs * kself.wk = list(self.w_ti) + list(self.w_si)fw = self.compute_fw()fw_list.append(fw)min_index = fw_list.index(min(fw_list))deta_wti = gt * lr * min_indexdeta_wsi = gs * lr * min_indexself.w_ti = temp_w_ti - deta_wtiself.w_si = temp_w_si - deta_wsiself.wk = list(self.w_ti) + list(self.w_si)def main():y = [[1,2,2],[2,1,1],[1,1,1],[1,2,2],[2,2,2],[2,1,2],[1,1,2],[1,2,1],[1,2,2],[1,2,2],[1,2,2],[1,2,2]]x = [[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1],[1,1,1]]CRF_test = CRF(y=y,y_num=2,N=3,x=x)print(CRF_test.Marix)print(CRF_test.compute_p(y[0]))print(CRF_test.predict_Py_i(2,1))print(CRF_test.Viterbi(x[0]))CRF_test.fit(50, how='IIS')print(CRF_test.wk)CRF_test.fit(500,how='GD',lr=0.001)print(CRF_test.wk)if __name__ == '__main__':main()#---------result-------------------- usr/bin/python3 /Users/zhengyanzhao/PycharmProjects/tongjixuexi/shixian2/CRF.py [[[2.71828183 1.64872127][1. 1. ]][[4.05519997 4.48168907][6.04964746 1.64872127]][[2.22554093 4.48168907][6.04964746 2.01375271]][[2.71828183 1. ][2.71828183 1. ]]] 0.06486783542907915 0.5082915274316868 (4.3, [1, 2, 1]) when iter is 11 shoulian [0.18764402126707205, 0.23235841506946808, 0.21344554371583913, 0.20112964244479178, 0.24547679607078648, 0.20600431503056774, 0.3669918383125717, 0.39593684451308797, 0.385216372799946, 0.377295537944856, 0.3971223364931556, 0.3658935638684794] when iter is 64 shoulian [0.23036225378918423, 0.0846981674992883, 0.0935981916969745, 0.03945814811020889, 0.0812481327750673, 0.16101306581281896, 0.43884913282239046, 0.15771695572586986, 0.37418378149743614, 0.2182414425846152, 0.16699382351353245, 0.4110185294594026]Process finished with exit code 0

總結

以上是生活随笔為你收集整理的统计学习方法第十一章作业：随机条件场—概率计算问题、IIS/GD学习算法、维特比预测算法 代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。