Python中如何为加权value-weight数据拟合分布并获取残差平方和？

一、手动实现加权拟合与残差计算 #

Scipy的scipy.stats系列拟合方法（比如norm.fit()）确实不直接返回误差指标，但完全可以手动实现加权拟合，再计算残差平方和，步骤如下：

1. 计算加权分布参数 #

• 正态分布：直接用加权均值和加权标准差作为分布参数

  import numpy as np
  from scipy.stats import norm, lognorm
  
  # 替换成你的真实数据
  values = np.array([11.12, 12.112])
  weights = np.array([0.001, 0.05])
  
  # 计算加权均值和标准差
  weighted_mean = np.average(values, weights=weights)
  weighted_var = np.average((values - weighted_mean)**2, weights=weights)
  weighted_std = np.sqrt(weighted_var)
  

• 对数正态分布：先对数值取对数，再计算加权参数

  log_values = np.log(values)
  log_weighted_mean = np.average(log_values, weights=weights)
  log_weighted_std = np.sqrt(np.average((log_values - log_weighted_mean)**2, weights=weights))
  

2. 计算加权残差平方和 #

残差平方和是拟合PDF与实际权重的差异加权平方和，注意要先把PDF缩放至和权重量级匹配（因为PDF积分是1，而权重总和也是1）：

# 正态分布残差平方和
norm_pdf = norm.pdf(values, loc=weighted_mean, scale=weighted_std)
# 用梯形法计算PDF积分，再缩放
pdf_integral = np.trapz(norm_pdf, values)
norm_pdf_scaled = norm_pdf / pdf_integral
weighted_sse_norm = np.sum(weights * (norm_pdf_scaled - weights)**2)

# 对数正态分布残差平方和
lognorm_pdf = lognorm.pdf(values, s=log_weighted_std, scale=np.exp(log_weighted_mean))
lognorm_integral = np.trapz(lognorm_pdf, values)
lognorm_pdf_scaled = lognorm_pdf / lognorm_integral
weighted_sse_lognorm = np.sum(weights * (lognorm_pdf_scaled - weights)**2)

二、用scipy.optimize.curve_fit做通用加权拟合 #

如果要适配更多分布，或者想让拟合更灵活，可以用curve_fit指定权重：

from scipy.optimize import curve_fit

# 定义正态分布PDF的拟合函数
def norm_pdf_func(x, mu, sigma):
    return norm.pdf(x, loc=mu, scale=sigma)

# 传入权重的倒数作为sigma参数，实现加权拟合
popt, _ = curve_fit(norm_pdf_func, values, weights, sigma=1/np.sqrt(weights), absolute_sigma=True)
fit_mu, fit_sigma = popt

# 计算残差平方和
fit_pdf = norm_pdf_func(values, fit_mu, fit_sigma)
fit_pdf_scaled = fit_pdf / np.trapz(fit_pdf, values)
weighted_sse = np.sum(weights * (fit_pdf_scaled - weights)**2)

三、关键注意点 #

• 权重是占比时，必须缩放PDF使其积分接近1，否则残差计算的量级会完全不对。
• 如果你的数据是分组区间的中点值，用np.trapz()计算PDF积分比简单的间距乘法更准确。
• 对数正态分布的scipy实现里，s参数对应对数转换后的标准差，scale是exp(对数均值)，别搞混参数。

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