统计学习方法第九章作业：三硬币EM算法、GMM高维高斯混合模型 代码实现

三硬币EM算法

import numpy as np
import math

class Three_coin:
    def __init__(self,pai=0.0,p=0.0,q=0.0):
        self.pai = pai
        self.p = p
        self.q = q

    def 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 = 0
        for x in X:
            result += math.log(self.comput_y_sita(x))
        return result

    def 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 np

class Gausian_EM:

    def __init__(self,Y,k):
        self.k = k
        self.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.001

    def 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/denominator

    def compute_log_likelihood(self):
        result = 0
        for 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 result

    def 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)/denominator
            for 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')
                break

    def 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]