基于直方图使用极大似然法求解混合指数分布的参数

任务内容 #

• 生成两个具有不同λ的随机指数分布，混合后用矩估计法提取两个λ
• 针对不同分箱数的直方图，通过极大似然法从混合数据集求解两个λ：
  • 对每个区间的概率密度函数（PDF）积分得到概率
  • 绘制极值附近的极大似然函数

当前卡点 #

无法基于直方图实现极大似然法求解两个λ，编写函数后未得到有效结果。

现有代码 #

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from random import choice
from scipy.optimize import fsolve
# Defining initial conditions for the random exponential distribution
n=100
lambda1=10
lambda2=50
# Creating two datasets with different lambdas
dataset1=pd.DataFrame(np.random.default_rng().exponential(scale=(1/lambda1),size=n))
dataset2=pd.DataFrame(np.random.default_rng().exponential(scale=(1/lambda2), size=n))
# Mixing all the data
mixeddata=pd.concat([dataset1,dataset2])
# Creating histograms for the different amount of bins
bins=[5,10,25,50]
fig, (ax1,ax2,ax3,ax4)= plt.subplots(4)
ax1.hist(mixeddata, bins=5,edgecolor='black', density=True)
ax2.hist(mixeddata, bins=10,edgecolor='black', density=True)
ax3.hist(mixeddata, bins=25,edgecolor='black', density=True)
ax4.hist(mixeddata, bins=50,edgecolor='black', density=True)
ax1.set_ylim(0,1)
ax2.set_ylim(0,1)
ax3.set_ylim(0,1)
ax4.set_ylim(0,1)

plt.show()
# Using method of moments to find lambda1 and lambda2 from the mixed dataset
sample_mean = mixeddata.mean().values[0]
sample_variance = mixeddata.var().values[0]

def equations(params):
    lambda1, lambda2, p = params
    eq1 = p * (1/lambda1) + (1-p) * (1/lambda2) - sample_mean
    eq2 = p * (1/(lambda1**2)) + (1-p) * (1/(lambda2**2)) + p*(1-p) * ((1/lambda1) - (1/lambda2))**2 - sample_variance
    # Assuming equal mixing probabilities for simplicity
    eq3 = p - 0.5
    return [eq1, eq2, eq3]

# Initial guesses for lambda1, lambda2, and p
initial_guesses = [10, 50, 0.5]
solutions = fsolve(equations, initial_guesses)

lambda1_est, lambda2_est, p_est = solutions
print(f'Estimated Parameters: lambda1 = {lambda1_est}, lambda2 = {lambda2_est}, p = {p_est}')

基于直方图的极大似然法实现方案 #

核心逻辑 #

直方图每个区间$[a_i, b_i)$的概率由混合指数分布积分得到：
$$P_i = p \cdot \int_{a_i}^{b_i} \lambda_1 e^{-\lambda_1 x} dx + (1-p) \cdot \int_{a_i}^{b_i} \lambda_2 e^{-\lambda_2 x} dx$$
似然函数为各区间观测频数$n_i$对应的多项分布似然，取对数后转为最小化负对数似然问题，通过数值优化求解参数。

实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import expon

# 生成混合数据
n = 100
lambda1_true = 10
lambda2_true = 50
rng = np.random.default_rng()
dataset1 = rng.exponential(scale=1/lambda1_true, size=n)
dataset2 = rng.exponential(scale=1/lambda2_true, size=n)
mixeddata = np.concatenate([dataset1, dataset2])

# 遍历不同分箱数
bin_counts = [5, 10, 25, 50]
for bins in bin_counts:
    # 计算直方图频数与区间边界
    counts, bin_edges = np.histogram(mixeddata, bins=bins, density=False)
    total_samples = len(mixeddata)

    # 定义负对数似然函数
    def neg_log_likelihood(params):
        lambda1, lambda2, p = params
        # 参数合法性检查
        if lambda1 <= 1e-3 or lambda2 <= 1e-3 or p <= 1e-3 or p >= 1-1e-3:
            return np.inf
        
        log_likelihood = 0.0
        for i in range(len(bin_edges)-1):
            a, b = bin_edges[i], bin_edges[i+1]
            # 计算单指数分布在区间内的概率
            prob1 = expon.cdf(b, scale=1/lambda1) - expon.cdf(a, scale=1/lambda1)
            prob2 = expon.cdf(b, scale=1/lambda2) - expon.cdf(a, scale=1/lambda2)
            # 混合分布概率
            mix_prob = p * prob1 + (1-p) * prob2
            # 避免对数为负无穷
            if mix_prob <= 1e-10:
                return np.inf
            log_likelihood += counts[i] * np.log(mix_prob)
        
        return -log_likelihood

    # 初始猜测（可使用矩估计结果或真实值）
    initial_guess = [lambda1_true, lambda2_true, 0.5]
    # 带边界约束的优化求解
    result = minimize(neg_log_likelihood, initial_guess, 
                     method='L-BFGS-B', 
                     bounds=[(1e-3, None), (1e-3, None), (1e-3, 1-1e-3)])
    
    lambda1_ml, lambda2_ml, p_ml = result.x
    print(f"分箱数 {bins}:")
    print(f"极大似然估计结果：lambda1={lambda1_ml:.2f}, lambda2={lambda2_ml:.2f}, p={p_ml:.2f}")
    print("-" * 50)

    # 绘制极值附近的负对数似然曲线（以lambda1为例）
    lambda1_range = np.linspace(lambda1_ml * 0.8, lambda1_ml * 1.2, 100)
    neg_ll_values = [neg_log_likelihood([l, lambda2_ml, p_ml]) for l in lambda1_range]

    plt.figure(figsize=(8, 4))
    plt.plot(lambda1_range, neg_ll_values, label='负对数似然')
    plt.scatter(lambda1_ml, neg_log_likelihood([lambda1_ml, lambda2_ml, p_ml]), 
               color='red', marker='o', label='极值点')
    plt.xlabel('lambda1')
    plt.ylabel('负对数似然值')
    plt.title(f'分箱数 {bins}：lambda1附近的负对数似然函数')
    plt.legend()
    plt.grid(alpha=0.3)
    plt.show()

关键说明 #

• 使用L-BFGS-B优化方法并设置参数边界，确保λ为正、混合比例p在(0,1)区间内
• 加入参数合法性检查与概率值下限判断，避免计算过程中出现无效值
• 每个分箱数对应独立的极大似然估计，并绘制极值附近的似然函数曲线，直观展示参数最优解的位置

内容的提问来源于stack exchange，提问作者dutchrunner