Gamma 函数拟合 OptimizeWarning：无法在 Python 中估计参数的协方差

Question

我正在尝试使用 scipy 中的 curve_fit 将函数拟合到直方图，但是运行程序时会出现警告：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.special import gamma as gamma_func
import tkinter as tk


def gaussian(x, m, s):
    return (1 / (s * np.sqrt(2 * np.pi))) * (np.exp((-1 / 2) * ((x - m) / s) ** 2))


def gamma(x, k, theta):
    return (x ** (k - 1) * np.exp(-x / theta)) / (theta ** k * gamma_func(k))


def display_gaussian(data_set, mean, sigma):
    plt.figure()
    y, bin_edges = np.histogram(data_set, bins=100, density=True)
    x = (bin_edges[:-1] + bin_edges[1:]) / 2
    plt.hist(data_set, density=True, color='skyblue', label='Histogram', edgecolor='black')

    init_gaussian = [mean, sigma]
    params_gaussian, cov_gaussian = curve_fit(gaussian, x, y, init_gaussian)


    fit = np.linspace(min(data_set), max(data_set), 1000)

    plt.plot(fit, gaussian(fit, *params_gaussian), 'r--',
             label=f'Fit: mean={params_gaussian[0]:.2f}, sigma={params_gaussian[1]:.2f}')
    plt.xlabel('Distribution')
    plt.ylabel('Frequency')
    plt.title('Gaussian distribution')
    plt.legend()
    plt.show(block=False)


def display_gamma(data_set, sigma):
    plt.figure()
    y, bin_edges = np.histogram(data_set, bins=100, density=True)
    x = (bin_edges[:-1] + bin_edges[1:]) / 2
    plt.hist(data_set, density=True, color='forestgreen', label='Histogram', edgecolor='black')

    shape_guess = 2.0
    init_gamma = [shape_guess, sigma]
    params_gamma, cov_gamma = curve_fit(gamma, x, y, init_gamma)


    fit = np.linspace(min(data_set), max(data_set), 1000)

    plt.plot(fit, gamma(fit, *params_gamma), 'y--',
             label=f'Fit: mean={params_gamma[0]:.2f}, sigma={params_gamma[1]:.2f}')
    plt.xlabel('Distribution')
    plt.ylabel('Frequency')
    plt.title('Gamma distribution')
    plt.legend()
    plt.show(block=False)


def run():
    mean = float(m_ent.get())
    sigma = float(s_ent.get())

    data_set = np.random.normal(mean, sigma, 1000)

    display_gaussian(data_set, mean, sigma)
    display_gamma(data_set, sigma)


root = tk.Tk()
root.title("Input parameters")
root.minsize(400, 300)

m_l = tk.Label(root, text="mean", font=15)
m_l.pack(pady=10)
m_ent = tk.Entry(root)
m_ent.pack(pady=10)

s_l = tk.Label(root, text="sigma", font=15)
s_l.pack(pady=10)
s_ent = tk.Entry(root)
s_ent.pack(pady=10)

run_btn = tk.Button(root, text="Run", command=run, font=15)
run_btn.pack(pady=20)

root.mainloop()

OptimizeWarning：无法估计参数的协方差

params_gamma, cov_gamma = curve_fit(gamma, x, y, init_gamma)

我不确定这里出了什么问题。

我认为问题出在直方图的箱上，并尝试将其设置为“自动”，同时更改密度，但警告再次出现。拟合高斯时也出现同样的问题，添加

density=True

后，成功估计了协方差，但它不适用于拟合伽玛。

Answer 1

您面临三个挑战：

• 将正态分布数据（可能为负）与具有正支持的伽马分布进行拟合；
• 在您有权访问原始数据时将 PDF 拟合到分箱数据（直方图），并且可以执行 MLE；
• 天真地编写 PDF，同时您可以依赖

  scipy.stats

  中实现的数字稳定版本。

这里是一个如何绘制伽玛样本并用 MLE 拟合它的示例：

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

law = stats.gamma(a=1.5, scale=0.8)
data = law.rvs(1_000)

def likelihood_factory(data):
    def wrapped(parameters):
        return - np.sum(stats.gamma(a=parameters[0], scale=parameters[1]).logpdf(data))
    return wrapped

likelihood = likelihood_factory(data)

solution = optimize.minimize(likelihood, x0=[1., 1.], tol=1e-4)
#  message: Optimization terminated successfully.
#  success: True
#   status: 0
#      fun: 1126.2309539053672
#        x: [ 1.508e+00  7.876e-01]
#      nit: 8
#      jac: [-1.526e-05 -3.052e-05]
# hess_inv: [[ 3.781e-03 -1.978e-03]
#            [-1.978e-03  1.447e-03]]
#     nfev: 33
#     njev: 11