scipy 曲线与自定义函数的直方图不匹配

n8ghc7c1  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(150)

我试图用一条曲线拟合一个直方图,但是得到的曲线是平坦的,尽管直方图不是平坦的。我该如何正确地拟合曲线?
我的当前代码:

import numpy as np

import matplotlib.pyplot as plt

from scipy.optimize import curve_fit

import pandas as pd

import scipy.optimize as optimization

x = np.random.normal(1e-10, 1e-7, size=10000)

def func(x, a, b, c):
    return a * (np.exp(-b*(x-c)**2))

bins=25

logbins = np.logspace(np.log10(1.0E-10),np.log10(1E-07),bins)

bin_heights, bin_borders, _ = plt.hist(x, bins=logbins, edgecolor='black', color='b')
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

x0    = np.array([0.0, 0.0, 0.0])

popt,cov = optimization.curve_fit(func, bin_centers, bin_heights,x0)
a,b,c=popt

popt, pcov = curve_fit(func, bin_centers, bin_heights, p0=[a,b,c])

x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 1000)

plt.plot(x_interval_for_fit, func(x_interval_for_fit, *popt), label='Fit',color='r')

plt.xscale('log')

plt.show()

结果是:

brccelvz

brccelvz1#

你得到的结果很差,因为你用来拟合直方图的函数看起来根本不像直方图的形状。通过使用一个简单的二阶插值函数,结果要好得多(尽管你可能会说不理想):

def func(x, a, b, c):
    return a * x**2 + b * x + c  # Simple 2nd-order polynomial

将它与您的代码一起使用(您可以删除两个优化步骤,并且只执行一次),我得到了以下结果:

在代码中,有一件事可能是无意的,那就是尽管创建了正态分布,但在直方图上,您决定以一种令人惊讶的方式对它们进行分组:只考虑分布的一端(由于您从1 e-10开始,并且您的分布居中在1 e-10中),并且您向右以对数方式增加条柱大小,您将在较大的条柱上得到更多的点。此外,您将忽略一半以上的点(那些小于1 e-10的,检查histdocumentation)。你可以通过以下方法检查:np.sum(bin_heights),您将看到计数小于len(x)的一半。
如果我上面所说的都是无意的,而你真的想找出随机数的原始生成高斯(但使用的是计数而不是密度),您应该直接在未修改的直方图上执行此操作。(在雅可比行列式具有广泛变化的值的意义上,取决于你去的维度)。你可以对“有多难”有一种感觉优化器通过手动计算误差并观察其变化情况来了解每个参数的影响由于您正在进行最小二乘优化,因此成本函数如下所示:

f = lambda a, b, c: sum((bin_heights - func(bin_centers, a, b, c))**2)

如果你在0的初始条件下使用这个函数,看看优化器看到了什么,你会注意到为什么收敛是困难的:

f(0, 0, 0)   # 8_407_566
f(0, 0, 1)   # 8_407_566, i.e. no change
f(0, 1, 0)   # 8_407_566, i.e. no change
f(1, 0, 0)   # 8_387_591, clearly an improvement

这不是一个数学上严格的论证(这来自于当ab接近于0时雅可比行列式的病态),但它给出了一个很好的直觉,说明哪里出了问题。
因此,你需要一个好的“第一猜测”。当然,在这个问题中,既然你知道你试图拟合一个高斯函数,你可以使用平均值、最大值和标准差作为初始参数(在你的例子中,是标准差平方的倒数)。下面的代码显示了一组“较少调整”的初始参数,它们也是有效的:

def func(x, a, b, c):  # Original function
    return a * np.exp(- b * (x - c)**2)

bins=25  # Let hist do its magic
bin_heights, bin_borders, _ = plt.hist(x, bins=bins, edgecolor='black', color='b')
bin_centers = bin_borders[:-1] + np.diff(bin_borders) / 2

x0 = np.array([600.0, 1.0e7, 0.0])

popt, cov = curve_fit(func, bin_centers, bin_heights, x0)

x_interval_for_fit = np.linspace(bin_borders[0], bin_borders[-1], 1000)

plt.plot(x_interval_for_fit, func(x_interval_for_fit, *x0), label='Initial guess',color='y')
plt.plot(x_interval_for_fit, func(x_interval_for_fit, *popt), label='Fit',color='r')
plt.show()

结果,包括由初始猜测创建的曲线(黄色)和插值的结果(红色):

相关问题