如何从Beta分布PDF数据拟合scipy.stats.beta的形状参数a与b？

问题背景 #

现有一段生成Beta分布概率密度函数（PDF）的代码，已知参数a=1.25、b=1.5，生成了对应x值的PDF输出数据。现在需要逆向操作：仅通过已知的x值和对应的PDF输出数据，拟合得到参数a和b。

生成PDF的代码如下：

import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1)

a = 1.25
b = 1.5

x = np.linspace(ss.beta.ppf(0.00, a, b),
                ss.beta.ppf(1, a, b), 100)

ax.plot(x, ss.beta.pdf(x, a, b),
       'r-', lw=5, alpha=0.6, label='beta pdf')

plt.show()

已知的PDF输出数据（对应上述x值）：

[0.         0.63154416 0.74719517 0.82263393 0.87936158 0.92490497
 0.96287514 0.9953117  1.02349065 1.04826858 1.07025112 1.08988379
 1.10750482 1.12337758 1.13771153 1.15067624 1.16241105 1.17303196
 1.18263662 1.19130806 1.19911747 1.20612637 1.21238825 1.21794995
 1.22285265 1.22713276 1.23082258 1.23395089 1.23654342 1.23862319
 1.24021091 1.24132519 1.2419828  1.2421989  1.24198713 1.24135982
 1.24032811 1.23890203 1.23709059 1.23490192 1.23234325 1.22942106
 1.22614106 1.2225083  1.21852712 1.21420131 1.209534   1.20452779
 1.19918472 1.19350628 1.18749347 1.18114672 1.17446601 1.16745076
 1.16009989 1.1524118  1.14438438 1.13601495 1.12730028 1.11823656
 1.10881937 1.09904366 1.0889037  1.07839304 1.06750447 1.05622997
 1.0445606  1.03248649 1.01999669 1.00707911 0.99372038 0.97990571
 0.96561876 0.95084139 0.93555349 0.91973269 0.90335407 0.88638972
 0.86880835 0.85057469 0.83164884 0.81198535 0.79153224 0.77022959
 0.74800782 0.7247854  0.70046587 0.67493372 0.64804878 0.61963821
 0.58948474 0.55730896 0.52274112 0.48527403 0.44417862 0.39833786
 0.34587496 0.2831383  0.20072303 0.        ]

尝试的拟合方法与错误结果 #

最初编写了仅使用两个点的拟合函数，随后修改为使用所有数据点的版本，但运行后得到完全错误的参数值：

a: 4.2467303147231366e-10
b: 1.9183210434237443

修改后的全数据点拟合代码如下：

def find_parameters(x, p):

    def objective(v):
        (a, b) = v
        return np.sum((ss.beta.cdf(x, a, b) - p) ** 2.0)

    # arbitrary initial guess of (1, 1) for parameters
    xopt = optimize.fmin(objective, (1, 1))
    return [xopt[0], xopt[1]]

fitted_a, fitted_b = find_parameters(x=x, p=ss.beta.pdf(x, a, b))

问题核心原因 #

代码中犯了根本性错误：将PDF数据当成了CDF数据来拟合。目标函数里用ss.beta.cdf()计算累积分布函数值，却和输入的PDF值做误差平方和，两者概念完全不同，导致优化器收敛到错误参数。

修正后的拟合代码 #

将目标函数中的ss.beta.cdf()替换为ss.beta.pdf()，同时优化导入逻辑：

import numpy as np
import scipy.stats as ss
from scipy.optimize import fmin

def find_parameters(x, p):
    def objective(v):
        a, b = v
        # 计算拟合的PDF值，和输入的p做误差平方和
        fitted_pdf = ss.beta.pdf(x, a, b)
        return np.sum((fitted_pdf - p) ** 2.0)

    # 初始猜测用(1,1)，也可根据PDF峰值位置调整更合理的初始值
    xopt = fmin(objective, (1, 1))
    return xopt[0], xopt[1]

# 假设x和p是已知的输入数据（即之前给出的x和PDF数组）
fitted_a, fitted_b = find_parameters(x=x, p=pdf_data)
print(f"拟合得到的a: {fitted_a}, b: {fitted_b}")

补充优化建议 #

1. 使用带约束的高效优化器：fmin是单纯形法优化器，可替换为scipy.optimize.minimize，选择支持边界约束的算法（Beta分布的a、b必须大于0），避免出现接近0的参数值：

from scipy.optimize import minimize

def find_parameters(x, p):
    def objective(v):
        a, b = v
        fitted_pdf = ss.beta.pdf(x, a, b)
        return np.sum((fitted_pdf - p) ** 2.0)

    # 设置参数边界：a>0，b>0
    bounds = ((1e-5, None), (1e-5, None))
    result = minimize(objective, (1, 1), method='L-BFGS-B', bounds=bounds)
    return result.x[0], result.x[1]

2. 利用Scipy内置拟合函数：Scipy提供ss.beta.fit()函数，可直接对样本数据拟合；若要拟合PDF曲线，手动编写目标函数更灵活，但内置函数可作为结果验证的参考。

验证结果 #

修正后的代码运行后，拟合参数会非常接近真实值a=1.25、b=1.5，例如可能得到：

a: 1.2498, b: 1.4997

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