多维高斯混合模型参数估计（python）

多维高斯混合分布EM算法参数估计（附python代码）

在博客高斯混合模型（R语言）的基础上，通过python实现了三维高斯混合模型的参数估计。

import numpy as np
import math
from scipy.linalg import *
from scipy import stats
import random
import sympy as sy


# Q函数计算
def Q_value(gama, alpha, X, phi):
    n= np.shape(X)[0]
    k = np.shape(alpha)[0]
    q= 0

    for j in range(n):
        for kk in range(k):
            q = q + alpha[kk][kk] * (np.log(alpha[kk][kk] + 1e-6) - 0.5 * np.log(phi[j][kk] + 1e-6))
    return q


def GM_EM(X, alpha, mu, List, ddmax, eposilon):
    # X          ：N*P维矩阵
    # alpha      ：K阶对角矩阵，对角线元素为a1,...,aK
    # mu         ：K*P阶矩阵，每一个行向量为P维向量μ1,...,μK
    # list       ：所有初始协方差矩阵
    # ddmax      ：最大迭代次数
    # eposilon   ：收敛阈值
    # plotnumber :是否绘图以及迭代第几次绘图
    # lastplot   :收敛时是否绘图

    K = np.shape(alpha)[1]
    N, P = np.shape(X)
    E = np.eye(P)
    a = (2 * math.pi) ** (-P / 2)  # a=2Π^(-P/2)
    Q = np.zeros((1, ddmax + 1))
    Q[0][0] = 1e-10
    cov = np.cov(X)

    i = 0
    while i < ddmax:
        # E步：
        # （1）求解条件期望E（γjk）构成的矩阵gama：维数为N*K
        # 1) 求phi：phi 为N*K维的矩阵，其第（i,j）元为N(xi|muj, sigmaj)

        # 构造向量c(|sigma1|, |sigma2|,...,|sigmaK|)
        det_sigma = np.zeros((1, K))
        for kk in range(K):
            det_sigma[0][kk] = det(List[kk] + 1e-6 * E)
        # 构造矩阵phi(N*K)
        phi = np.zeros((N, K))
        for j in range(N):
            for kk in range(K):
                x = (X[j] - mu[kk]).reshape((1, P))
                phi[j][kk] = stats.multivariate_normal(mu[kk], List[kk]).pdf(X[j])

        # 2) 构造γ_hat矩阵，维数为N*K
        # gama_pre维数为N*K，其第(i, j)元为αj*N(xi,muj,sigmaj)
        gama_pre = np.dot(phi, alpha)
        gama = np.zeros((N,K))
        # gama矩阵即为gama_pre的每一行除以其行和，矩阵维数为N*K
        gama_rowsum = np.sum(gama_pre, 1)
        for ii in range(N):
            gama[ii] = gama_pre[ii] / gama_rowsum[ii]

        # 记录Q函数值
        Q[0][i+1] = Q_value(gama, alpha, X, phi)

        # M步
        # 各参数极大似然估计

        # gama矩阵第k列的列和
        gama_colsum  = np.sum(gama, 0)
        # print(gama_colsum)
        # (1) alpha是K阶对角矩阵，对角线元素为α1,...,αK，其中αk等于gama矩阵的每一列的列和再除以N
        alpha = np.diag(gama_colsum/N)
        print(alpha)
        # (2) mu是K*P阶矩阵，每一个行向量为P维向量μ1,...,μK
        for ii in range(K):
            mu[ii] = np.dot(gama.T, X)[ii]/gama_colsum[ii]

        # (3) list包括所有协方差矩阵，有Σ1,...,ΣK组成
        for kk in range(K):
            # mu_single为P*N的矩阵，矩阵每一列都是相同的均值向量μk
            mu_single = []
            for ii in range(N):
                mu_single.append(mu[kk])
            mu_single = np.array(mu_single).T
            x1 = np.diag(gama[:,kk])
            List[kk] = np.dot(np.dot((X.T-mu_single),np.diag(gama[:,kk])),(X.T-mu_single).T)/gama_colsum[kk]
        # 制图
        # A = phi/a
        # P_pre= np.dot(A,alpha)
        # # P为预测的概率矩阵
        # P = P_pre/np.sum(P_pre, 0)
        #
        # variance =
        minus = np.abs(Q[0][i+1]-Q[0][i])
        print(minus)

        if minus < eposilon:
            break
        else:
            i += 1
        print(i)
    return alpha, mu, List


N = 500
K = 2
P = 3
a1 = 0.1
mu1 = [2, 2, 2]
sigma1 = np.diag([1, 1, 1])

a2 = 0.9
mu2 = [1.2, 1.2, 1.2]
sigma2 = np.diag([1.44, 1.44, 1.44])

X = np.zeros((N, P))
u = np.random.random(size=N)
fenbu = 0

for i in range(N):
    if u[i] < a1:
        X[i] = np.random.multivariate_normal(mu1, sigma1)
        fenbu += 1
    if u[i] >= a1:
        X[i] = np.random.multivariate_normal(mu2, sigma2)

alpha = np.diag([0.2,0.8])
mu = [[1.2,1, 1],[1.4,1.2, 1.2]]
mu = np.array(mu)
List = []
for i in range(K):
    if i == 0:
        List.append(np.diag([1, 1, 1]))
    else:
        List.append((np.diag([1.3,1.3, 1.3])))

List = np.array(List)
ddmax = 3000

alpha_n, mu_n, List_n = GM_EM(X, alpha, mu,List, ddmax, eposilon=1e-6)
print('alpha = ', alpha_n)
print('mu = ', mu_n[1])
print('list = ', List_n[1])