基于Scipy.stats为已知概率值拟合分布及解决curve_fit参数问题

已知概率0.1对应值1，0.5对应值4，0.9对应值7，拟合正态分布的代码示例：

正态分布的中位数等于均值，因此直接取中间值4作为均值。再利用分位数与标准正态分位数的关系计算标准差：

import scipy.stats as st

# 已知分位数与对应概率
quantiles = [1, 4, 7]
probs = [0.1, 0.5, 0.9]

# 正态分布参数计算
mu = quantiles[1]  # 中位数即均值
z_01 = st.norm.ppf(probs[0])  # 标准正态0.1分位数
z_09 = st.norm.ppf(probs[2])  # 标准正态0.9分位数
sigma = (quantiles[2] - quantiles[0]) / (z_09 - z_01)

print(f"拟合得到的正态分布参数：均值μ={mu:.4f}，标准差σ={sigma:.4f}")

# 验证拟合效果
for q, p in zip(quantiles, probs):
    pred_p = st.norm.cdf(q, loc=mu, scale=sigma)
    print(f"值{q}的预测概率：{pred_p:.4f}，目标概率：{p}")

报错核心原因是curve_fit无法自动识别不同分布的参数数量，且未提供合理的初始参数猜测，部分分布参数存在约束条件（如必须为正）。以下是修正后的代码，包含异常捕获、动态参数处理、RMSE评估：

import numpy as np
import scipy.optimize as opt
import scipy.stats as st

# 你的数据：数值对应累积概率（可替换为你原代码中的[0.114, 0.28, 0.32]）
data = np.array([1, 4, 7])
prob = np.array([0.1, 0.5, 0.9])

# 待测试的分布列表（移除trapz，它是积分函数不是概率分布）
DISTRIBUTIONS = [
    st.alpha, st.anglit, st.arcsine, st.beta, st.betaprime, st.bradford, st.burr, st.burr12, st.cauchy, st.chi, st.chi2, st.cosine,
    st.erlang, st.expon, st.exponnorm, st.exponweib, st.exponpow, st.f, st.fatiguelife, st.fisk,
    st.foldnorm, st.genlogistic, st.genexpon, st.genpareto,
    st.genextreme, st.gausshyper, st.gamma, st.gengamma, st.genhalflogistic, st.gilbrat, st.gompertz, st.gumbel_r,
    st.gumbel_l, st.halfcauchy, st.halflogistic, st.halfnorm, st.halfgennorm, st.hypsecant, st.invgamma, st.invgauss,
    st.invweibull, st.johnsonsb, st.johnsonsu, st.ksone, st.kstwobign, st.laplace, st.levy, st.levy_l,
    st.logistic, st.loggamma, st.loglaplace, st.lognorm, st.lomax, st.maxwell, st.mielke, st.nakagami, st.ncx2, st.ncf,
    st.nct, st.norm, st.pareto, st.pearson3, st.powerlaw, st.powerlognorm, st.powernorm, st.rdist, st.reciprocal, st.norminvgauss,
    st.rayleigh, st.rice, st.recipinvgauss, st.semicircular, st.t, st.triang, st.truncexpon, st.truncnorm, st.tukeylambda,
    st.uniform, st.vonmises, st.vonmises_line, st.wald, st.weibull_min, st.weibull_max, st.wrapcauchy
]

# 存储最优拟合结果
best_fit = None
best_rmse = float('inf')

for dist in DISTRIBUTIONS:
    try:
        # 获取分布的形状参数数量，加上loc和scale为总参数数
        num_shape_params = len(dist.shapes.split(',')) if dist.shapes else 0
        total_params = num_shape_params + 2

        # 设置初始参数猜测：形状参数设为1，loc用中位数，scale用数据范围的1/6
        initial_guess = [1.0]*num_shape_params + [np.median(data), np.ptp(data)/6]

        # 定义拟合用的cdf函数，拆分参数为形状参数、loc、scale
        def cdf_func(x, *params):
            shape_params = params[:num_shape_params]
            loc, scale = params[-2], params[-1]
            scale = np.abs(scale)  # 确保尺度参数为正
            return dist(*shape_params, loc=loc, scale=scale).cdf(x)

        # 设置参数约束：形状参数非负，scale最小为1e-6避免为0
        bounds = ([0.0]*num_shape_params + [-np.inf, 1e-6], 
                  [np.inf]*num_shape_params + [np.inf, np.inf])
        params, _ = opt.curve_fit(cdf_func, data, prob, p0=initial_guess, bounds=bounds)

        # 计算RMSE评估拟合效果
        pred_probs = cdf_func(data, *params)
        rmse = np.sqrt(np.mean((pred_probs - prob)**2))

        print(f"分布{dist.name}：参数={params.round(4)}，RMSE={rmse:.6f}")

        # 更新最优拟合记录
        if rmse < best_rmse:
            best_rmse = rmse
            best_fit = (dist.name, params, rmse)

    except Exception as e:
        # 跳过拟合失败的分布
        print(f"分布{dist.name}拟合失败：{str(e)}")
        continue

print("\n最优拟合结果：")
print(f"分布：{best_fit[0]}，参数={best_fit[1].round(4)}，RMSE={best_fit[2]:.6f}")

关键修正点说明： #

• 移除无效项：删除trapz，它是积分工具而非概率分布
• 动态参数适配：根据不同分布的形状参数数量，自动构建初始猜测和拟合函数
• 参数约束：对尺度参数设置下界，避免出现非正值导致分布无效
• 异常处理：跳过参数不收敛、约束冲突等拟合失败的情况
• 量化评估：用RMSE衡量预测概率与目标概率的偏差，筛选最优分布

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