三参数Weibull分布拟合双列数据的问题及初始值求解咨询

问题 #

现有数据第一列为x值（孔径），第二列为y值（PSD），需要用三参数Weibull分布拟合。尝试scipy.stats.exponweib和scipy.stats.weibull_min均失败，自行定义函数用scipy.optimize.curve_fit拟合，但调整初始值后仍无法得到满意结果。代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Load data
poresize, psd, psd_std = np.loadtxt("data.txt", unpack=True)

# Define the Weibull distribution function with the requested form
def weibull_func(x, a, b, c):
    return a * b * (x - c) ** (b - 1) * np.exp(-a * (x - c) ** b)

# Perform curve fitting
popt, pcov = curve_fit(weibull_func, poresize, psd, p0=[1, 1, 1])

# Plot the original data and the fitted curve
plt.scatter(poresize, psd, label='Data')
x_range = np.linspace(min(poresize), max(poresize), 100)
plt.plot(x_range, weibull_func(x_range, *popt), 'r-', label='Fitted curve (Weibull)')
plt.xlabel('Particle Size')
plt.ylabel('PSD')
plt.title('Fitting Weibull Distribution')
plt.legend()
plt.grid(True)
plt.show()

# Display the optimized parameters
a_opt, b_opt, c_opt = popt
print("Optimized Parameters:")
print("a:", a_opt)
print("b:", b_opt)
print("c:", c_opt)

函数形式改编自Matlab的Weibull分布（将x替换为x-c），请问代码是否有根本性错误？如何获取合适的初始猜测值？

分析与解决 #

1. 代码的根本性问题 #

• 定义域未处理：三参数Weibull要求x > c（位置参数），当x <= c时，(x - c)为非正数，若b-1不是整数会产生复数或NaN，直接导致拟合崩溃。需在函数中添加判断，当x <= c时返回0。
• 参数定义与标准形式不匹配：你的函数是a*b*(x-c)^(b-1)*exp(-a*(x-c)^b)，而标准三参数Weibull PDF的形式为：
  $$f(x) = \frac{b}{a^b} (x-c)^{b-1} \exp\left(-\left(\frac{x-c}{a}\right)^b\right)$$
  两者参数对应关系不一致，会让拟合算法难以收敛，也导致初始值的猜测逻辑混乱。

2. 快速获取初始值的方法 #

方法一：利用scipy.stats.weibull_min先拟合转换

scipy.stats.weibull_min原生支持三参数Weibull，先通过它得到初始参数，再转换为你定义的函数参数：

• weibull_min.fit返回的参数为(形状参数b, 位置参数c, 尺度参数a)
• 转换关系：你的函数中a = 1/(scale^b)，b = 形状参数，c = 位置参数

方法二：手动估算

• 位置参数c：c必须小于所有x值，可取数据中最小x值的90%（如c = np.min(poresize)*0.9）
• 形状参数b：对数据做对数变换，先计算CDF（累加PSD并归一化），再对ln(-ln(1-CDF))和ln(x-c)做线性拟合，斜率即为b的初始值
• 尺度相关参数a：用b和c的初始值代入，通过最小二乘法粗略估算a

3. 修改后的示例代码 #

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

# Load data
poresize, psd, psd_std = np.loadtxt("data.txt", unpack=True)

# 修正后的Weibull函数，处理定义域问题
def weibull_func(x, a, b, c):
    # 当x <= c时返回0，避免复数/NaN
    mask = x > c
    result = np.zeros_like(x)
    x_valid = x[mask]
    result[mask] = a * b * (x_valid - c) ** (b - 1) * np.exp(-a * (x_valid - c) ** b)
    return result

# Step1: 用weibull_min获取初始参数
# weibull_min.fit返回 (形状参数, 位置参数, 尺度参数)
shape_init, loc_init, scale_init = weibull_min.fit(poresize, floc=np.min(poresize)*0.8)
# 转换为自定义函数的初始参数
a_init = 1 / (scale_init ** shape_init)
b_init = shape_init
c_init = loc_init

# Step2: 用curve_fit拟合
popt, pcov = curve_fit(weibull_func, poresize, psd, p0=[a_init, b_init, c_init])

# 绘图
plt.scatter(poresize, psd, label='Data')
x_range = np.linspace(min(poresize), max(poresize), 100)
plt.plot(x_range, weibull_func(x_range, *popt), 'r-', label='Fitted curve (Weibull)')
plt.xlabel('Particle Size')
plt.ylabel('PSD')
plt.title('Fitting Weibull Distribution')
plt.legend()
plt.grid(True)
plt.show()

# 输出参数
a_opt, b_opt, c_opt = popt
print("Optimized Parameters:")
print("a:", a_opt)
print("b:", b_opt)
print("c:", c_opt)

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