Hiểu mô hình hỗn hợp Gaussian

Question

Tôi đang cố gắng hiểu kết quả từ việc thực hiện mô hình hỗn hợp gaussian scikit-learn. Hãy xem ví dụ sau:

#!/opt/local/bin/python 
import numpy as np 
import matplotlib.pyplot as plt 
from sklearn.mixture import GaussianMixture 

# Define simple gaussian 
def gauss_function(x, amp, x0, sigma): 
    return amp * np.exp(-(x - x0) ** 2./(2. * sigma ** 2.)) 

# Generate sample from three gaussian distributions 
samples = np.random.normal(-0.5, 0.2, 2000) 
samples = np.append(samples, np.random.normal(-0.1, 0.07, 5000)) 
samples = np.append(samples, np.random.normal(0.2, 0.13, 10000)) 

# Fit GMM 
gmm = GaussianMixture(n_components=3, covariance_type="full", tol=0.001) 
gmm = gmm.fit(X=np.expand_dims(samples, 1)) 

# Evaluate GMM 
gmm_x = np.linspace(-2, 1.5, 5000) 
gmm_y = np.exp(gmm.score_samples(gmm_x.reshape(-1, 1))) 

# Construct function manually as sum of gaussians 
gmm_y_sum = np.full_like(gmm_x, fill_value=0, dtype=np.float32) 
for m, c, w in zip(gmm.means_.ravel(), gmm.covariances_.ravel(), 
       gmm.weights_.ravel()): 
    gmm_y_sum += gauss_function(x=gmm_x, amp=w, x0=m, sigma=np.sqrt(c)) 

# Normalize so that integral is 1  
gmm_y_sum /= np.trapz(gmm_y_sum, gmm_x) 

# Make regular histogram 
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=[8, 5]) 
ax.hist(samples, bins=50, normed=True, alpha=0.5, color="#0070FF") 
ax.plot(gmm_x, gmm_y, color="crimson", lw=4, label="GMM") 
ax.plot(gmm_x, gmm_y_sum, color="black", lw=4, label="Gauss_sum") 

# Annotate diagram 
ax.set_ylabel("Probability density") 
ax.set_xlabel("Arbitrary units") 

# Draw legend 
plt.legend() 
plt.show() 

đây đầu tiên tôi tạo ra một phân phối mẫu xây dựng từ Gaussian, sau đó phù hợp với một mô hình hỗn hợp gaussian đến những dữ liệu này. Tiếp theo, tôi muốn tính toán xác suất cho một số đầu vào đã cho. Thuận tiện, việc triển khai scikit cung cấp phương thức score_samples để thực hiện điều đó. Bây giờ tôi đang cố gắng hiểu những kết quả này. Tôi luôn nghĩ rằng, tôi chỉ có thể lấy các thông số của gaussians từ GMM phù hợp và xây dựng phân phối rất giống nhau bằng cách tổng hợp chúng và sau đó chuẩn hóa tích phân thành 1. Tuy nhiên, như bạn thấy trong cốt truyện, các mẫu được rút ra từ phương pháp score_samples phù hợp hoàn hảo (đường màu đỏ) với dữ liệu gốc (biểu đồ màu xanh), phân phối được tạo thủ công (đường màu đen) thì không. Tôi muốn hiểu nơi mà suy nghĩ của tôi đã đi sai và tại sao tôi không thể tự mình xây dựng bản phân phối bằng cách tổng hợp các gaussians như được đưa ra bởi sự phù hợp của GMM!?! Cảm ơn rất nhiều cho bất kỳ đầu vào!

Answer 1

Chỉ trong trường hợp bất cứ ai trong tương lai đang băn khoăn về những điều tương tự: Người ta phải bình thường hóa các thành phần riêng lẻ, không phải là số tiền:

import numpy as np 
import matplotlib.pyplot as plt 
from sklearn.mixture import GaussianMixture 

# Define simple gaussian 
def gauss_function(x, amp, x0, sigma): 
    return amp * np.exp(-(x - x0) ** 2./(2. * sigma ** 2.)) 

# Generate sample from three gaussian distributions 
samples = np.random.normal(-0.5, 0.2, 2000) 
samples = np.append(samples, np.random.normal(-0.1, 0.07, 5000)) 
samples = np.append(samples, np.random.normal(0.2, 0.13, 10000)) 

# Fit GMM 
gmm = GaussianMixture(n_components=3, covariance_type="full", tol=0.001) 
gmm = gmm.fit(X=np.expand_dims(samples, 1)) 

# Evaluate GMM 
gmm_x = np.linspace(-2, 1.5, 5000) 
gmm_y = np.exp(gmm.score_samples(gmm_x.reshape(-1, 1))) 

# Construct function manually as sum of gaussians 
gmm_y_sum = np.full_like(gmm_x, fill_value=0, dtype=np.float32) 
for m, c, w in zip(gmm.means_.ravel(), gmm.covariances_.ravel(), gmm.weights_.ravel()): 
    gauss = gauss_function(x=gmm_x, amp=1, x0=m, sigma=np.sqrt(c)) 
    gmm_y_sum += gauss/np.trapz(gauss, gmm_x) * w 

# Make regular histogram 
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=[8, 5]) 
ax.hist(samples, bins=50, normed=True, alpha=0.5, color="#0070FF") 
ax.plot(gmm_x, gmm_y, color="crimson", lw=4, label="GMM") 
ax.plot(gmm_x, gmm_y_sum, color="black", lw=4, label="Gauss_sum", linestyle="dashed") 

# Annotate diagram 
ax.set_ylabel("Probability density") 
ax.set_xlabel("Arbitrary units") 

# Make legend 
plt.legend() 

plt.show()