混合Gamma分布5个参数精确推断：解决curve_fit未拟合出双峰问题

已知离散数据点服从如下混合Gamma分布：

y = w * gamma.pdf(x, alpha1, scale=scale1) + (1-w) * gamma.pdf(x, alpha2, scale=scale2)

需要精确推断这5个参数（w, alpha1, scale1, alpha2, scale2）。使用curve_fit拟合时未得到双峰形态，希望能精准逼近参数，神经网络等方法也可接受。

用户提供的原始代码：

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import gaussian_kde,gamma
data1 = np.random.gamma(4, 1, 200)
data2 = np.random.gamma(6,2, 200)
data = np.concatenate((data1, data2))
x=np.linspace(0,np.max(data)+2,800)
y=gaussian_kde(data)(x) 
initial_params = [0.5, 2, 1, 2, 1]  
params, params_covariance = curve_fit(two_gamma, x, y, p0=initial_params, maxfev=50000)
w, alpha1, scale1, alpha2, scala2 = params
plt.figure(figsize=(10, 6))
sns.histplot(data, bins=20, kde=False, color='y', label='Data density', alpha=0.5, stat='probability')
plt.plot(x, y, marker='o', linestyle='', markersize=1, label='Data distribution')
y_fit=w*gamma.pdf(x, alpha1, scale=scale1)+(1-w)*gamma.pdf(x, alpha2, scale=scala2)
plt.plot(x, y_fit, 'r-', linewidth=1, alpha=0.7,label='Mixture gamma distribution')
plt.legend(fontsize=8, loc='upper right')
plt.title("Expression distribution of gamma mixture")
plt.xlabel("Expression")

问题修正与解决方案 #

1. 修复curve_fit的使用问题 #

原始代码存在三个核心问题：未定义拟合函数two_gamma、初始参数偏离真实值、用KDE结果拟合不如直接用原始数据直方图准确。修正后的代码如下：

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import gamma

# 生成数据
data1 = np.random.gamma(4, 1, 200)
data2 = np.random.gamma(6, 2, 200)
data = np.concatenate((data1, data2))

# 定义混合Gamma PDF拟合函数
def two_gamma(x, w, alpha1, scale1, alpha2, scale2):
    return w * gamma.pdf(x, alpha1, scale=scale1) + (1 - w) * gamma.pdf(x, alpha2, scale=scale2)

# 用原始数据的概率密度直方图作为拟合目标
hist, bins = np.histogram(data, bins=30, density=True)
x_hist = (bins[:-1] + bins[1:]) / 2  # 取区间中点

# 设置接近真实值的初始参数，添加参数约束
initial_params = [0.5, 3, 1, 5, 2]  
params, params_covariance = curve_fit(
    two_gamma, 
    x_hist, 
    hist, 
    p0=initial_params,
    bounds=([0, 0, 0, 0, 0], [1, np.inf, np.inf, np.inf, np.inf]),  # w∈[0,1]，其余参数>0
    maxfev=100000
)

w, alpha1, scale1, alpha2, scale2 = params
print(f"拟合参数：w={w:.2f}, alpha1={alpha1:.2f}, scale1={scale1:.2f}, alpha2={alpha2:.2f}, scale2={scale2:.2f}")

# 绘图展示
x = np.linspace(0, np.max(data)+2, 800)
y_fit = two_gamma(x, w, alpha1, scale1, alpha2, scale2)

plt.figure(figsize=(10, 6))
sns.histplot(data, bins=30, kde=False, color='y', alpha=0.5, stat='density', label='数据直方图')
plt.plot(x, y_fit, 'r-', linewidth=1.5, alpha=0.7, label='拟合混合Gamma分布')
plt.legend(fontsize=10, loc='upper right')
plt.title("混合Gamma分布拟合结果")
plt.xlabel("数值")
plt.show()

2. 更可靠的方法：EM算法（期望最大化） #

混合分布的参数推断，EM算法是标准解法，能有效避免局部最优问题，稳定性远超curve_fit。手动实现的简化版代码如下：

import numpy as np
from scipy.special import psi, polygamma
from scipy.stats import gamma

def em_mixture_gamma(data, n_components=2, max_iter=1000, tol=1e-6):
    n = len(data)
    # 初始化参数
    w = np.ones(n_components) / n_components
    alpha = np.random.uniform(2, 8, n_components)
    scale = np.random.uniform(1, 3, n_components)
    
    log_likelihood_prev = -np.inf
    for _ in range(max_iter):
        # E步：计算每个数据点属于两个分量的后验概率（责任度）
        gamma_pdfs = np.array([gamma.pdf(data, a=alpha[i], scale=scale[i]) for i in range(n_components)])
        responsibilities = (w * gamma_pdfs) / (np.sum(w * gamma_pdfs, axis=0) + 1e-10)
        
        # M步：更新参数
        # 更新混合权重
        w = np.sum(responsibilities, axis=1) / n
        
        # 用牛顿法更新每个Gamma分量的alpha和scale
        for i in range(n_components):
            resp = responsibilities[i]
            sum_resp = np.sum(resp)
            sum_resp_x = np.sum(resp * data)
            sum_resp_log_x = np.sum(resp * np.log(data))
            
            alpha_i = alpha[i]
            # 牛顿法迭代求解alpha的MLE
            for __ in range(10):
                psi_alpha = psi(alpha_i)
                polygamma_alpha = polygamma(1, alpha_i)
                scale_i = sum_resp_x / (sum_resp * alpha_i)
                log_scale = np.log(scale_i)
                numerator = sum_resp * (psi_alpha - log_scale) + sum_resp_log_x - sum_resp * np.log(scale_i)
                denominator = sum_resp * polygamma_alpha
                delta = numerator / denominator
                alpha_i += delta
                if abs(delta) < 1e-8:
                    break
            alpha[i] = alpha_i
            scale[i] = scale_i
        
        # 检查对数似然收敛
        log_likelihood = np.sum(np.log(np.sum(w * gamma_pdfs, axis=0) + 1e-10))
        if abs(log_likelihood - log_likelihood_prev) < tol:
            break
        log_likelihood_prev = log_likelihood
    
    return w, alpha, scale

# 使用EM算法拟合
w_em, alpha_em, scale_em = em_mixture_gamma(data)
print(f"EM算法拟合参数：w={w_em[0]:.2f}, alpha1={alpha_em[0]:.2f}, scale1={scale_em[0]:.2f}, alpha2={alpha_em[1]:.2f}, scale2={scale_em[1]:.2f}")

3. 神经网络方法（可选） #

针对大数据量或复杂分布场景，可使用神经网络拟合混合分布，通过负对数似然作为损失函数训练：

import torch
import torch.nn as nn
import torch.optim as optim
from scipy.stats import gamma

class GammaMixtureNet(nn.Module):
    def __init__(self):
        super().__init__()
        self.fc = nn.Sequential(
            nn.Linear(1, 64),
            nn.ReLU(),
            nn.Linear(64, 32),
            nn.ReLU(),
            nn.Linear(32, 5)
        )
    
    def forward(self, x):
        params = self.fc(x)
        w = torch.sigmoid(params[:, 0])  # 约束w∈[0,1]
        alpha1 = torch.exp(params[:, 1])  # 约束alpha>0
        scale1 = torch.exp(params[:, 2])
        alpha2 = torch.exp(params[:, 3])
        scale2 = torch.exp(params[:, 4])
        return w, alpha1, scale1, alpha2, scale2

# 准备数据
data_tensor = torch.tensor(data.reshape(-1, 1), dtype=torch.float32)

# 初始化模型与优化器
model = GammaMixtureNet()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# 训练过程
for epoch in range(5000):
    optimizer.zero_grad()
    w, a1, s1, a2, s2 = model(data_tensor)
    
    # 计算负对数似然损失
    pdf1 = torch.tensor(gamma.pdf(data_tensor.numpy(), a=a1.detach().numpy(), scale=s1.detach().numpy()), dtype=torch.float32)
    pdf2 = torch.tensor(gamma.pdf(data_tensor.numpy(), a=a2.detach().numpy(), scale=s2.detach().numpy()), dtype=torch.float32)
    mixture_pdf = w * pdf1 + (1 - w) * pdf2
    loss = -torch.mean(torch.log(mixture_pdf + 1e-10))
    
    loss.backward()
    optimizer.step()
    
    if epoch % 500 == 0:
        print(f"Epoch {epoch}, Loss: {loss.item():.4f}")

# 获取拟合参数
w_nn, a1_nn, s1_nn, a2_nn, s2_nn = model(data_tensor)
print(f"神经网络拟合参数：w={w_nn.mean().item():.2f}, alpha1={a1_nn.mean().item():.2f}, scale1={s1_nn.mean().item():.2f}, alpha2={a2_nn.mean().item():.2f}, scale2={s2_nn.mean().item():.2f}")

内容的提问来源于stack exchange，提问作者NiguLas张三