三参数Weibull分布拟合曲线与数据存间隙及拟合效果差的问题排查

问题描述 #

将数据拟合至三参数Weibull分布时遇到两个问题：

• 原始数据直方图与拟合PDF曲线存在明显间隙
• 使用初始参数[1,1,1]时拟合效果极差

代码片段 #

import numpy as np  
from scipy.optimize import minimize  
from scipy.special import gamma  
import matplotlib.pyplot as plt  
from scipy.stats import weibull_min

def weibull_log_likelihood(params, data):
    shape, scale, loc = params
    log_likelihood = -np.sum(weibull_min.logpdf(data, shape, loc=loc, scale=scale))
    return -log_likelihood  

def estimate_weibull_params(data):
    initial_params = [shape, scale, loc] 
    bounds=[(1, None), (1, None), (1, None)]
    result = minimize(weibull_log_likelihood, initial_params, args=(data,), method='nelder-mead', bounds=bounds) 
    return result.x  

def weibull_pdf(x, shape, scale, loc):
    return (shape / scale) * ((x - loc) / scale) ** (shape - 1) * np.exp(-((x - loc) / scale) ** shape)

shape = 7.5
scale = 150
loc = 350
size = 100
data = weibull_min.rvs(shape, loc=loc, scale=scale, size=size)
   
estimated_params = estimate_weibull_params(data)
shape, scale, loc = estimated_params

print(f"Estimated Parameters: Shape = {shape}, Scale = {scale}, Location = {loc}")

x = np.arange(1000)
pdf = weibull_pdf(x, shape, scale, loc)  
plt.hist(data, bins=20, density=True, alpha=0.6, color='g')  
plt.plot(x, pdf, 'r-', lw=2)  
plt.xlim(0, 1000)  
plt.show()  

错误分析与修正方案 #

1. 对数似然函数逻辑完全颠倒 #

当前代码中，对数似然函数的计算逻辑错误：

log_likelihood = -np.sum(weibull_min.logpdf(data, shape, loc=loc, scale=scale))
return -log_likelihood

这相当于返回原始对数似然的和，而我们需要最大化对数似然，使用minimize时应该返回负的对数似然（最小化负对数似然等价于最大化对数似然）。正确写法应为：

def weibull_log_likelihood(params, data):
    shape, scale, loc = params
    return -np.sum(weibull_min.logpdf(data, shape, loc=loc, scale=scale))

2. 初始参数依赖全局变量+边界设置不合理 #

• estimate_weibull_params函数中initial_params = [shape, scale, loc]依赖全局变量，代码耦合性高；若替换为[1,1,1]，与真实参数（shape=7.5, scale=150, loc=350）差距过大，Nelder-Mead无梯度优化易陷入局部最优。
• loc的边界设为(1, None)错误，三参数Weibull的支撑集是x >= loc，loc不能超过数据最小值，否则对数似然会出现负无穷。

修正方案：

def estimate_weibull_params(data, initial_params=None):
    if initial_params is None:
        # 用数据特征初始化，更接近真实值
        loc_init = np.min(data) - 1
        scale_init = np.std(data)
        shape_init = 2
        initial_params = [shape_init, scale_init, loc_init]
    # loc上限设为数据最小值，保证所有数据在支撑集内
    bounds=[(1e-3, None), (1e-3, None), (-np.inf, np.min(data))]
    # 换用L-BFGS-B方法，边界处理更友好、收敛更稳定
    result = minimize(weibull_log_likelihood, initial_params, args=(data,), method='L-BFGS-B', bounds=bounds) 
    return result.x

3. 自定义PDF未处理无效区间 #

三参数Weibull在x <= loc时PDF值为0，但当前函数中x < loc时，(x-loc)为负数，计算幂次会得到复数或nan，导致曲线在有效区间外出现异常值，视觉上形成“间隙”。

修正方案：

def weibull_pdf(x, shape, scale, loc):
    mask = x > loc
    pdf = np.zeros_like(x, dtype=np.float64)
    x_valid = x[mask] - loc
    pdf[mask] = (shape / scale) * (x_valid / scale) ** (shape - 1) * np.exp(-(x_valid / scale) ** shape)
    return pdf

修正后完整代码 #

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

def weibull_log_likelihood(params, data):
    shape, scale, loc = params
    return -np.sum(weibull_min.logpdf(data, shape, loc=loc, scale=scale))

def estimate_weibull_params(data, initial_params=None):
    if initial_params is None:
        loc_init = np.min(data) - 1
        scale_init = np.std(data)
        shape_init = 2
        initial_params = [shape_init, scale_init, loc_init]
    bounds=[(1e-3, None), (1e-3, None), (-np.inf, np.min(data))]
    result = minimize(weibull_log_likelihood, initial_params, args=(data,), method='L-BFGS-B', bounds=bounds)
    return result.x

def weibull_pdf(x, shape, scale, loc):
    mask = x > loc
    pdf = np.zeros_like(x, dtype=np.float64)
    x_valid = x[mask] - loc
    pdf[mask] = (shape / scale) * (x_valid / scale) ** (shape - 1) * np.exp(-(x_valid / scale) ** shape)
    return pdf

# 生成模拟数据
true_shape = 7.5
true_scale = 150
true_loc = 350
size = 100
data = weibull_min.rvs(true_shape, loc=true_loc, scale=true_scale, size=size)

# 用[1,1,1]作为初始参数测试
estimated_params = estimate_weibull_params(data, initial_params=[1,1,1])
shape, scale, loc = estimated_params

print(f"真实参数: Shape = {true_shape}, Scale = {true_scale}, Location = {true_loc}")
print(f"估计参数: Shape = {shape:.2f}, Scale = {scale:.2f}, Location = {loc:.2f}")

# 绘制图表
x = np.linspace(0, 1000, 1000)
pdf = weibull_pdf(x, shape, scale, loc)
plt.hist(data, bins=20, density=True, alpha=0.6, color='g', label='原始数据直方图')
plt.plot(x, pdf, 'r-', lw=2, label='拟合PDF曲线')
plt.xlim(true_loc - 50, true_loc + 3*true_scale)  # 聚焦数据有效区间
plt.legend()
plt.show()

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