EM算法的python實現

本文參考自：https://www.jianshu.com/p/154ee3354b59 和 李航博士的《統計學習方法》

1.

 

2. 創建觀測結果數據 

def createData(m,n):
    y = np.mat(np.zeros((m,n)))
    
    for i in range(m):
        for j in range(n):
            # 通過產生隨機數，每一行表示一次實驗結果 
            y[i,j] = random.randint(0,1)
    return y

輸出一下，觀察一下結果：

data = createData(1,10)  #一共做了三次實驗，每次觀測到了10個硬幣C出現的結果
data

　結果： matrix([[0., 0., 1., 1., 1., 1., 0., 1., 0., 1.]]) 

3.  EM算法的實現過程

 

def EM(arr_y,theta,tol,num_iter):
    #初始化參數
    PI = 0 
    P = 0 
    Q = 0 
    m,n = np.shape(arr_y)
    mat_y = arr_y.getA()  #返回的是一個numpy array 的數組
    
    for i in range(num_iter):
        miu = []
        PI = np.copy(theta[0])  # 深拷貝
        P = np.copy(theta[1])
        Q = np.copy(theta[2])
        for j in range(m):
            miu_value = (PI*(P**mat_y[j]) *((1-P)**(1-mat_y[j]))) / \
            (PI*(P**mat_y[j])*((1-P)**(1-mat_y[j])) + (1-PI)*(Q**mat_y[j])*((1-Q)**(1-mat_y[j])))
            miu.append(miu_value)
            
        sum1 = 0.0 
        for j in range(m):
            sum1 += miu[j]
        theta[0] = sum1 / m 
        
        sum1 = 0.0 
        sum2 = 0.0 
        for j in range(m):
            sum1 += miu[j] * mat_y[j]
            sum2 += miu[j]
        theta[1] = sum1 / sum2
        
        sum1 = 0.0 
        sum2 = 0.0 
        for j in range(m):
            sum1 += (1-miu[j])* mat_y[j]
            sum2 += (1-miu[j])
        theta[2] = sum1 / sum2
        
        print("-----------------------------")
        print(theta)
        if (abs(theta[0] - PI) <= tol and abs(theta[1] - P) <= tol 
            and abs(theta[2] - Q <= tol)):
            print("迭代完畢，參數已經收斂")
            break 
    return PI,P,Q 

4. 主函數的實現 （注意：這里的輸入數據（與《統計學習方法》的輸入數據一樣））

if __name__ == "__main__":
    mat_y = np.mat(np.zeros((10, 1)))
    mat_y[0,0] = 1
    mat_y[1,0] = 1
    mat_y[2,0] = 0
    mat_y[3,0] = 1
    mat_y[4,0] = 0
    mat_y[5,0] = 0
    mat_y[6,0] = 1
    mat_y[7,0] = 0
    mat_y[8,0] = 1
    mat_y[9,0] = 1
    theta = [0.5, 0.5, 0.5]
    print(mat_y)
    PI,P,Q = EM(mat_y,theta,0.001,100)
    print(PI,P,Q)

#本文的輸出結果

[[1.]
 [1.]
 [0.]
 [1.]
 [0.]
 [0.]
 [1.]
 [0.]
 [1.]
 [1.]]
-----------------------------
[array([0.5]), array([0.6]), array([0.6])]
-----------------------------
[array([0.5]), array([0.6]), array([0.6])]
迭代完畢，參數已經收斂
[0.5] [0.6] [0.6] 

和書上的輸出結果是一樣的