使用Scipy拟合Rice分布出现异常结果的问题排查

你的代码存在几个关键错误，导致拟合结果无意义，以下是问题点和修正方案：

1. 自定义Rice PDF公式核心缺失 #

你实现的rice_pdf遗漏了零阶第一类修正贝塞尔函数I₀，这是Rice分布PDF的核心组成部分。正确的Rice分布PDF公式为：

$f(x; b, \nu) = \frac{x}{b^2} \exp\left(-\frac{x^2 + \nu^2}{2b2}\right) I_0\left(\frac{x\nu}{b^2}\right)$

同时你添加的amplitude参数是冗余的——直方图已经设置density=True，PDF本身是归一化的，不需要额外振幅缩放。

2. 拟合后分布对象的参数逻辑错误 #

你创建拟合分布时使用rice(nu, loc=scale, scale=np.sqrt(b**2 + scale**2))完全不符合scipy.stats.rice的参数定义，该分布的构造格式为rice(b, loc=loc, scale=scale)，其中b是形状参数（对应$\nu/\sigma$），loc是位置偏移，scale是缩放参数$\sigma$。

3. 初始猜测参数不合理 #

初始猜测的参数数量冗余且取值不符合Rice分布的参数逻辑，容易导致拟合算法收敛到局部最优或无意义解。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import rice
from scipy.optimize import curve_fit
from scipy.special import iv  # 导入零阶修正贝塞尔函数

# 符合scipy定义的Rice PDF函数
def rice_pdf(x, b, scale, loc=0):
    x_shifted = x - loc
    # 处理Rice分布定义域（x >= loc）
    x_shifted[x_shifted < 0] = 0
    nu = b * scale
    return (x_shifted / scale**2) * np.exp(-(x_shifted**2 + nu**2) / (2 * scale**2)) * iv(0, (x_shifted * nu) / scale**2)

def fit_rice_distribution_to_histogram(hist_data, bins):
    bin_centers = (bins[:-1] + bins[1:]) / 2

    # 基于数据特征的合理初始猜测
    data_std = np.std(bin_centers)
    initial_guess = [0.8, data_std, 0]
    # 限制参数范围避免无意义解
    bounds = ([0, 0, -np.inf], [np.inf, np.inf, min(bin_centers)])
    
    params, covariance = curve_fit(rice_pdf, bin_centers, hist_data, p0=initial_guess, bounds=bounds)
    b_fit, scale_fit, loc_fit = params

    # 创建正确的拟合分布对象
    fitted_distribution = rice(b=b_fit, loc=loc_fit, scale=scale_fit)
    return fitted_distribution, b_fit, scale_fit, loc_fit

if __name__ == "__main__":
    # 真实分布参数
    nu_true = 8.5
    sigma_true = 10.5
    sample_size = 1000  # 增大样本量让直方图更接近真实PDF
    b_true = nu_true / sigma_true

    # 生成Rice分布数据
    data = rice.rvs(b=b_true, scale=sigma_true, size=sample_size)

    # 绘制归一化直方图
    hist_data, bins, _ = plt.hist(data, bins=30, density=True, alpha=0.5, label="生成数据")
    plt.xlabel("数值")
    plt.ylabel("概率密度")

    # 拟合分布
    fitted_dist, b_fit, scale_fit, loc_fit = fit_rice_distribution_to_histogram(hist_data, bins)

    # 绘制拟合曲线
    x = np.linspace(min(bins), max(bins), 1000)
    pdf_values = fitted_dist.pdf(x)
    plt.plot(x, pdf_values, 'r-', label="拟合Rice分布")
    plt.legend()
    plt.show()

    # 输出参数对比
    print(f"真实参数：b={b_true:.4f}, sigma={sigma_true}, nu={nu_true}")
    print(f"拟合参数：b={b_fit:.4f}, scale={scale_fit:.4f}, loc={loc_fit:.4f}")
    print(f"推导的nu（b*scale）：{b_fit*scale_fit:.4f}")

额外优化说明 #

• 增大样本量到1000，让直方图更接近真实PDF，提升拟合稳定性；
• 添加参数边界限制，避免拟合出负数等无意义参数；
• 移除冗余的amplitude参数，保持PDF与归一化直方图的一致性；
• 修正了分布对象的参数传递逻辑，确保拟合曲线与真实分布匹配。

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