高斯混合模型(GMM)和EM演算法 —— python實現

小··明發表於2024-03-27
模型演算法Python

一、EM演算法

EM演算法是一種迭代演算法，用於含有隱含變數的機率模型引數的極大似然估計。設Y為觀測隨機變數的資料，Z為隱藏的隨機變數資料，Y和Z一起稱為完全資料。

觀測資料的似然函式為：
image

模型引數θ的極大似然估計為：
image

這個問題只有透過迭代求解，下面給出EM演算法的迭代求解過程：

step1、選擇合適的引數初值θ(0)，開始迭代

step2、E步，求期望。θ(i)為第i次迭代θ的估計值，在第i+1步，計算下面的Q函式：
image

Q函式為logP(Y,Z|θ)關於在給定觀測資料Y和當前引數θ(i)下對隱藏變數Z的條件機率分佈P(Z|Y,θ(i))的期望。

step3、M步，求極大化。求使Q函式極大化的θ，確定第i+1次迭代的引數估計：

step4、重複第2、3步，直到收斂。

EM演算法對初值的選取比較敏感，且不能保證找到全域性最優解。

二、在高斯混合模型（GMM）中的應用

一維高斯混合模型：
image

多維高斯混合模型：
image

wk（k=1，2，……，K）為混合項係數，和為1。∑為協方差矩陣。θ=（wk，uk，σk）。

設有N個可觀測資料yi，它們是這樣產生的：先根據機率wk選擇第k個高斯分佈模型，生成觀測資料yi。yi是已知的，但yi屬於第j個模型是未知的，是隱藏變數。用Zij表示隱藏變數，含義是第i個資料屬於第j個模型的機率。先寫出完全資料的似然函式，然後確定Q函式，要最大化期望，對wk、uk、σk求偏導並使其為0。可得高斯混合模型引數估計的EM演算法（以高維資料為例）：

step1、引數賦初始值，開始迭代

step2、E步，計算混合項係數Zij的期望 E[Zij]：
image

step3、M步，計算新一輪迭代的引數模型：
image

image

image

step4、重複第2、3步，直到收斂。

前文搬運自連結：https://blog.csdn.net/u014157632/article/details/65442165

三、python程式示例

首先初始化分佈。

MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];

並用這兩個分佈各生成1000個散點。

點選檢視程式碼 
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 初始化分佈引數
MU1 = np.array([1, 2])
SIGMA1 = np.array([[1, 0], [0, 0.5]])
MU2 = np.array([-1, -1])
SIGMA2 = np.array([[1, 0], [0, 1]])

# 生成資料點
data1 = np.random.multivariate_normal(MU1, SIGMA1, 1000)
data2 = np.random.multivariate_normal(MU2, SIGMA2, 1000)

X = np.concatenate([data1, data2])  # [2000, 2]

#對X進行隨機打亂，此步復現時不可忽略
np.random.shuffle(X)

接下來繪製二維散點圖。

點選檢視程式碼 
plt.scatter(data1[:, 0], data1[:, 1], label='Class 1')
plt.scatter(data2[:, 0], data2[:, 1], label='Class 2')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scatter Plot of Data Points')
plt.legend()
plt.show()

image

EM演算法實現高斯混合模型的擬合

點選檢視程式碼 

import math


def fit_gaussian_mixture_model(X, n_components=2, max_iter=100):
    # 合併資料
    gmm = GaussianMixture(n_components=n_components, max_iter=max_iter)
    gmm.fit(X)

    return gmm.means_, gmm.covariances_, gmm.weights_


###擬合函式！
def my_fit_GMM(X, n_components=2, max_iter=1000):

    # 給資料排序
    n_samples, n_features = X.shape
    X_SortIndex = np.lexsort((X[:, 0], X[:, 1]))
    newArr = []
    for id in X_SortIndex:
        newArr.append(X[id])
    X = np.array(newArr)

    # 初始均值向量與協方差矩陣
    mu = X[:n_components, :]                  # [n_components, n_features]
    # cov = np.array(([[1.,0.],[0.,1.]],[[1.,0.],[0.,1.]]))
    cov = np.ones(shape=(n_components, n_features, n_features))        # [n_components, n_features, n_features]
    interval = n_samples // n_components

    for i in range(n_components):
        cov[i, ...] = np.cov(np.transpose(X[i*interval: (i+1)*interval, ...]))
        mu[i, ...] = np.mean(X[i*interval: (i+1)*interval, ...], axis=0)
    weights = np.ones(shape=(n_components, )) / n_components                      # [n_components, ]


    # 直接編寫即可，不一定要使用結構體。
    expectation = np.zeros(shape=(n_samples, n_components))
    old_mu = np.copy(mu)

    X = np.matrix(X)
    mu = np.matrix(mu)
    expectation = np.matrix(expectation)

    for iteration in range(max_iter):
        old_mu = np.copy(mu)
        # 更新期望E
        for i in range(n_samples):
            for j in range(n_components):
                expectation[i, j] = math.exp(-(X[i,:]-mu[j,:])*np.linalg.inv(cov[j,...])*np.transpose(X[i,:]-mu[j,:])/2)\
                                    / np.sqrt(np.linalg.det(cov[j,...]))
                expectation[i, j] *= weights[j]

        t = np.sum(expectation, axis=1).reshape(-1, 1)
        expectation = expectation / t

        # 更新模型（均值向量）
        sum_E = np.zeros(shape=(n_components, ))
        for j in range(n_components):
            sum_E[j] = np.sum(expectation[:, j])
            mu[j, ...] = np.sum(np.multiply(expectation[:, j], X), axis=0) / sum_E[j]

        # 更新模型（協方差矩陣），混合項係數
        for j in range(n_components):
            cov_s = np.matrix(np.zeros((n_features, n_features)))
            for i in range(n_samples):
                cov_i = (X[i, ...] - mu[j, ...]).T * (X[i, ...] - mu[j, ...])
                cov_s = cov_s + np.multiply(expectation[i, j], cov_i)
            cov[j, ...] = np.array(cov_s) / sum_E[j]
            weights[j] = sum_E[j] / n_samples


        # 停止迭代的條件
        if np.abs(np.linalg.det(old_mu - mu)) < 1e-6:
            break

    means = np.array(mu)
    return means, cov, weights

# 請改為用my_fit_GMM擬合散點分佈。並輸出估計的均值,方差和混合項係數。

# means, cov, weights = fit_gaussian_mixture_model(X)
#
means, cov, weights = my_fit_GMM(X)

print("Data1.Means:")
print(means[0])
print("Data2.Means:")
print(means[1])

print("Data1.Covariances:")
print(cov[0])
print("Data2.Covariances:")
print(cov[1])

print("Weights:")
print(weights)

Data1.Means:
[-0.98994072 -0.99780233]
Data2.Means:
[1.00491397 1.98848228]
Data1.Covariances:
[[ 1.00762612 -0.04130926]
[-0.04130926 0.93607269]]
Data2.Covariances:
[[1.07037303 0.04512588]
[0.04512588 0.47164631]]
Weights:
[0.49222433 0.50777567]

點選檢視程式碼 
#繪製結果的機率密度函式，用於驗證演算法效果

from scipy.stats import multivariate_normal

def draw_results(m,s,w):

    # 定義兩個二維高斯分佈的引數
    m1 = m[0]
    s1 = s[0]

    m2 = m[1]
    s2 = s[1]

    # 生成網格點
    x, y = np.mgrid[-5:5:.01, -5:5:.01]
    pos = np.dstack((x, y))

    # 計算每個點的機率密度值
    rv1 = multivariate_normal(m1, s1)
    rv2 = multivariate_normal(m2, s2)
    z1 = rv1.pdf(pos)
    z2 = rv2.pdf(pos)

    # 混合兩個高斯分佈的機率密度函式
    z = w[0] * z1 + w[1] * z2  # 設定混合比例

    # 繪製3D圖
    fig = plt.figure(figsize=(5, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, cmap='viridis')

    # 設定座標軸標籤和標題
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Probability Density')
    ax.set_title('3D Plot of Mixture Gaussian Probability Density')

    plt.show()

draw_results(means, cov, weights)
#看到擬合的視覺化結果了

image

相關文章