三硬幣EM算法

import numpy as np import mathclass Three_coin:def __init__(self,pai=0.0,p=0.0,q=0.0):self.pai = paiself.p = pself.q = qdef comput_y_sita(self,y):return self.pai*self.p**y*(1-self.p)**(1-y) + (1-self.pai)*self.q**y*(1-self.q)**(1-y)def log_P_Y_sita(self,X):result = 0for x in X:result += math.log(self.comput_y_sita(x))return resultdef compute_ui(self,y):return self.pai*self.p**y*(1-self.p)**(1-y)/self.comput_y_sita(y)def fit(self,X,max_iter):self.n = len(X)for i in range(max_iter):p_u1 = np.array([self.compute_ui(x)*x for x in X])ui = np.array([self.compute_ui(x) for x in X])q_ui = np.array([(1 - self.compute_ui(x))*x for x in X])self.pai = 1/self.n*sum(ui)self.p = sum(p_u1)/sum(ui)self.q = sum(q_ui)/sum(1-ui)def main():X = [1,1,0,1,0,0,1,0,1,1]Three_coin_ = Three_coin(0.46,0.55,0.67)Three_coin_.fit(X,10)print(Three_coin_.pai,Three_coin_.p,Three_coin_.q)if __name__ == '__main__':main() ###result######################## /usr/bin/python3 /Users/zhengyanzhao/PycharmProjects/tongjixuexi/shixian2/three_coin_model.py 0.461862835113919 0.5345950037850112 0.6561346417857326

GMM高維高斯混合模型

import math import numpy as npclass Gausian_EM:def __init__(self,Y,k):self.k = kself.Y = np.array(Y)self.feature_num = len(Y[0])self.N = len(Y)self.uk = []self.sik = []for i in range(k):self.uk.append(np.random.rand(self.feature_num))self.sik.append(np.random.rand(self.feature_num,self.feature_num))self.ak = np.array([1/k]*k)self.rjk = np.zeros((self.N,k)) + 0.001def caculate_y_sita(self,y,k_index):covdet = np.linalg.det(self.sik[k_index] + np.eye(self.feature_num) * 0.001)covinv = np.linalg.inv(self.sik[k_index] + np.eye(self.feature_num) * 0.001)denominator = ((2*math.pi)**self.feature_num * np.abs(covdet))**(1/2)numerator = np.exp(-0.5*((y-self.uk[k_index]).dot(covinv).dot(y-self.uk[k_index])))return numerator/denominatordef compute_log_likelihood(self):result = 0for y in self.Y:result += np.log(np.array(np.sum([self.caculate_y_sita(y,k)*self.ak[k] for k in range(self.k)])))return resultdef fit(self,max_iter):for iter in range(max_iter):log_likelihood = self.compute_log_likelihood()for n in range(self.N):denominator = np.sum([self.ak[k] * self.caculate_y_sita(self.Y[n], k) for k in range(self.k)])for k_index in range(self.k):self.rjk[n][k_index] = self.ak[k_index]*self.caculate_y_sita(self.Y[n], k_index)/denominatorfor k in range(self.k):self.ak[k] = np.sum([self.rjk[j][k] for j in range(self.N)])/float(self.N)self.sik[k] = np.sum([self.rjk[j][k]*((self.Y[j]-self.uk[k]).reshape(self.feature_num,1).dot((self.Y[j]-self.uk[k]).reshape(1,self.feature_num))) \for j in range(self.N)],axis=0)/np.sum([self.rjk[j][k] for j in range(self.N)])print([self.Y[j]*self.rjk[j][k] for j in range(self.N)])self.uk[k] = np.sum([self.Y[j]*self.rjk[j][k] for j in range(self.N)],axis=0)/np.sum([self.rjk[j][k] for j in range(self.N)])## print('---------------------------')# print(self.ak)# print(self.sik)# print(self.uk)new_log_likelihood = self.compute_log_likelihood()if new_log_likelihood - log_likelihood < 0.0001:print('small fit')breakdef predict(self,x):return np.argmax([self.caculate_y_sita(x,k) for k in range(self.k)])def main():X = []Y = []with open('../data/iris.data', 'r') as f:for i in f:data = i.split(',')X.append([float(j) for j in data[:4]])Y.append(data[4])Y = [1 if i == 'Iris-setosa\n' else 0 for i in Y]Gausian_EM_ = Gausian_EM(X,2)Gausian_EM_.fit(50)print([Gausian_EM_.predict(x) for x in X])print(Y)if __name__ == '__main__':main()#####result################################### …… --------------------------- [0.33332862 0.66667138] [array([[0.12176211, 0.09828544, 0.01581506, 0.01033654],[0.09828544, 0.14226048, 0.01144559, 0.01120903],[0.01581506, 0.01144559, 0.02950403, 0.00558421],[0.01033654, 0.01120903, 0.00558421, 0.01126411]]), array([[0.43497483, 0.12094219, 0.44886966, 0.16550395],[0.12094219, 0.10961752, 0.14138142, 0.07923294],[0.44886966, 0.14138142, 0.674851 , 0.28587703],[0.16550395, 0.07923294, 0.28587703, 0.17863612]])] [array([5.00600713, 3.41801571, 1.46400229, 0.24399922]), array([6.26198756, 2.871996 , 4.90597454, 1.67599028])] small fit [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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, 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, 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, 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, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

總結(jié)

以上是生活随笔為你收集整理的统计学习方法第九章作业：三硬币EM算法、GMM高维高斯混合模型 代码实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。