機器學習(五):混合高斯聚類GMM(求聚類標籤)+PCA降維(3維降2維)習題

孤飛發表於2023-04-08
  1. 使用混合高斯模型 GMM,計算如下資料點的聚類過程:
    \(Data = np.array([1,2,6,7])\)
    均值初值為:
    \(\mu_1, \mu_2 = 1, 5\)
    權重初值為:
    \(w_1, w_2 = 0.5, 0.5\)
    方差:
    \(std_1, std_2 = 1, 1\)
    \(K = 2\)
    10 次迭代後資料的聚類標籤是多少?

採用python程式碼實現:

from scipy import stats
import numpy as np

#初始化資料
Data = np.array([1,2,6,7])
w1 , w2 = 0.5, 0.5
mu1 , mu2 = 1, 5
std1 , std2 = 1, 1

n = len(Data) # 樣本長度
zij=np.zeros([n,2])
for t in range(10):
    # E-step 依據當前引數,計算每個資料點屬於每個子分佈的機率
    z1_up = w1 * stats.norm(mu1 ,std1).pdf(Data)
    z2_up = w2*stats.norm(mu2 , std2).pdf(Data)
    z_all = (w1*stats.norm(mu1 ,std1).pdf(Data)+w2*stats.norm(mu2 ,std2).pdf(Data))+0.001
    rz1 = z1_up/z_all # 為甲分佈的機率
    rz2 = z2_up/z_all # 為乙分佈的機率
    # M-step 依據 E-step 的結果,更新每個子分佈的引數。
    mu1 = np.sum(rz1*Data)/np.sum(rz1)
    mu2 = np.sum(rz2*Data)/np.sum(rz2)
    std1 = np.sum(rz1*np.square(Data-mu1))/np.sum(rz1)
    std2 = np.sum(rz2*np.square(Data-mu2))/np.sum(rz2)
    w1 = np.sum(rz1)/n
    w2 = np.sum(rz2)/n
for i in range(n):
    zij[i][0] = rz1[i]/(rz1[i]+rz2[i])
    zij[i][1] = rz2[i]/(rz1[i]+rz2[i])

labels = np.argmax(zij, axis=1)#輸出每一行的最大值,0或1  axis表示返回每一行中最大值所在列的索引
print(labels)

聚類標籤輸出結果:[0 0 1 1]

也就是說,10 次迭代後資料的聚類標籤是1,2歸為0類6,7歸為1

附註:

如果 axis 為 None,那麼 np.argmax 會將陣列展平為一維,然後返回最大值的索引。例如:

>>> a = np.array([[1, 2], [3, 4]])
>>> np.argmax(a)
3

如果 axis 為 0,那麼 np.argmax 會沿著第一個維度(行)進行最大值的查詢,返回每一列中最大值所在的行索引。例如:

>>> a = np.array([[1, 2], [3, 4]])
>>> np.argmax(a, axis=0)
array([1, 1])

如果 axis 為 1,那麼 np.argmax 會沿著第二個維度(列)進行最大值的查詢,返回每一行中最大值所在的列索引。例如:

>>> a = np.array([[1, 2], [3, 4]])
>>> np.argmax(a, axis=1)
array([1, 1])

在之前問題中,np.argmax([gamma1, gamma2], axis=0) 的意思是沿著第一個維度(gamma1 和 gamma2)進行最大值的查詢,返回每個資料點屬於哪個子分佈的機率更大。

  1. 假設我們的資料集有 10 個 3 維資料, 需要用 PCA 降到 2 維特徵。

    array([
        [ 3.25, 1.85, -1.29],
        [ 3.06, 1.25, -0.18],
        [ 3.46, 2.68, 0.64],
        [ 0.3 , -0.1 , -0.79],
        [ 0.83, -0.21, -0.88],
        [ 1.82, 0.99, 0.16],
        [ 2.78, 1.75, 0.51],
        [ 2.08, 1.5 , -1.06],
        [ 2.62, 1.23, 0.04],
        [ 0.83, -0.69, -0.61]])
    

    給出求解過程

解:

  1. 對所有的樣本進行中心化:

\[x^{(i)}=x^{(i)}-\frac{1}{m} \sum_{j=1}^{m} x^{(j)} \]

得到:

X=np.array([
     [ 1.147  0.825 -0.944]
     [ 0.957  0.225  0.166]
     [ 1.357  1.655  0.986]
     [-1.803 -1.125 -0.444]
     [-1.273 -1.235 -0.534]
     [-0.283 -0.035  0.506]
     [ 0.677  0.725  0.856]
     [-0.023  0.475 -0.714]
     [ 0.517  0.205  0.386]
     [-1.273 -1.715 -0.264]])
  1. 計算樣本的協方差矩陣 $X X^{T} $
covM2=np.array([[1.26344556 1.08743889 0.32030889], 
[1.08743889 1.11076111 0.31611111],
[0.32030889 0.31611111 0.45449333]])
  1. 對矩陣 $X X^{T} $ 進行特徵值分解

取出最大的 \(\mathrm{n}^{\prime}\) 個特徵值對應的特徵向量 $ \left(w_{1}, \ldots, w_{n^{\prime}}\right) $, 將所有的特徵向量標準化後,組成特徵向量矩陣 \(W\)

3.1求出特徵值:

eigval=np.array([2.38219729 0.09637041 0.35013229])

3.2特徵向量標準化:

eigvec=np.array([
[ 0.71144     0.67380165 -0.19961077],
[ 0.66498574 -0.73733944 -0.11884665],
[ 0.22725997  0.04818606  0.97264126]])

3.3取出特徵值最大的2個特徵值索引,也就是\([2.38, 0.35]\)對應的第1列和第3列:

indexes=[2 0]

3.4特徵向量矩陣W:(對eigvec取了第3列和第1列)

W=np.array([
[-0.19961077  0.71144   ], 
[-0.11884665   0.66498574], 
[ 0.97264126   0.22725997]])
  1. 對樣本集中的每一個樣本 \(x^{(i)}\) , 轉化為新的樣本 \(z^{(i)}=W^{T} x^{(i)}\) ,得到輸出樣本集 $D=\left(z^{(1)}, \ldots z^{(m)}\right) $

X:3×10 W:3×2 \(x\cdot W =10\times3 \quad 3\times2\) 因為輸入行列轉置,結果是一致的

D=np.array([
     [-1.24517539  1.15010151]
     [-0.05630956  0.86819503]
     [ 0.49146125  2.29005381]
     [ 0.06174799 -2.1317387 ]
     [-0.1185103  -1.84827733]
     [ 0.55280596 -0.10961848]
     [ 0.6112806   1.15829407]
     [-0.74632697  0.13724149]
     [ 0.24787719  0.5918589 ]
     [ 0.20114923 -2.10611029]])

程式碼:

import numpy as np

X=np.array([
    [ 3.25, 1.85, -1.29],
    [ 3.06, 1.25, -0.18],
    [ 3.46, 2.68, 0.64],
    [ 0.3 , -0.1 , -0.79],
    [ 0.83, -0.21, -0.88],
    [ 1.82, 0.99, 0.16],
    [ 2.78, 1.75, 0.51],
    [ 2.08, 1.5 , -1.06],
    [ 2.62, 1.23, 0.04],
    [ 0.83, -0.69, -0.61]])

def pca(X, d):
    # Centralization中心化
    means = np.mean(X, 0)
    X = X - means
    print(X)
    # Covariance Matrix 計算樣本協方差矩陣
    M=len(X)
    X=np.mat(X)    
    covM2=np.cov(X.T)
    # 求出特徵值,特徵值分解
    eigval , eigvec = np.linalg.eig(covM2)
    indexes = np.argsort(eigval)[-d:]
    W = eigvec[:, indexes]
    return X*W
print(pca(X, 2))

附註:

np.cov()是一個用於計算協方差矩陣的函式,它可以接受一個或兩個陣列作為引數,返回一個二維陣列,表示協方差矩陣。

協方差矩陣是一個對稱矩陣,它的對角線元素表示各個變數的方差,非對角線元素表示兩個變數之間的協方差。協方差反映了兩個變數的線性相關程度,如果協方差為正,說明兩個變數正相關;如果協方差為負,說明兩個變數負相關;如果協方差為零,說明兩個變數無相關性。

np.cov()的用法如下:

np.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None)

引數說明:

  • m: 一個一維或二維的陣列,表示多個變數和觀測值。如果是一維陣列,表示一個變數的觀測值;如果是二維陣列,每一行表示一個變數,每一列表示一個觀測值。
  • y: 可選引數,另一個一維或二維的陣列,表示另一組變數和觀測值,必須和m具有相同的形狀。
  • rowvar: 可選引數,布林值,預設為True。如果為True,表示每一行代表一個變數;如果為False,表示每一列代表一個變數。
  • bias: 可選引數,布林值,預設為False。如果為False,表示計算無偏協方差(除以n-1);如果為True,表示計算有偏協方差(除以n)。
  • ddof: 可選引數,整數,預設為None。如果不為None,則覆蓋由bias隱含的預設值。ddof=0表示計算有偏協方差;ddof=1表示計算無偏協方差。
  • fweights: 可選引數,一維陣列或整數,預設為None。表示每次觀測的頻率權重。
  • aweights: 可選引數,一維陣列,預設為None。表示每個變數的可靠性權重。

返回值:

  • 一個二維陣列,表示協方差矩陣。

舉例說明:

import numpy as np

# 生成兩組隨機資料
x = np.random.randn(10)
y = np.random.randn(10)

# 計算x和y的協方差矩陣
cov_xy = np.cov(x,y)
print(cov_xy)
# 輸出:
[[ 0.8136679  -0.01594772]
 [-0.01594772  0.84955963]]

# 計算x和y的相關係數矩陣
corr_xy = np.corrcoef(x,y)
print(corr_xy)
# 輸出:
[[ 1.         -0.01904402]
 [-0.01904402  1.        ]]

相關文章