我使用以下方法生成了随机数据：
bkg= 240-140*np.random.power(3.5,50000)
我通过使用将点绘制成直方图
h_all = plt.hist(all,bins=binedges,histtype='step')
我的问题是，如果我知道 pdf（在本例中称为“bkg”），我是否可以使用 scipy.optimize 生成一条与生成的点完美匹配的曲线，以及曲线的方程是什么？
【问题讨论】：
标签： python optimization
首先，请注意您的bkg 不是概率密度函数 (pdf)。相反，它是来自 pdf 的 观察 列表。通过在这个观察列表上调用matplotlib.pyplot.hist，您可以看到一条近似于概率密度函数（偏移和缩放版本）的曲线。如果给定了这条曲线，则可以很好地估计建模所需的参数，前提是您已经先验地获得了参数化模型。
例如：
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
offset, scale, a, nsamples = 240, -140, 3.5, 500000
bkg = offset + scale*np.random.power(a, nsamples) # values range between (offset, offset+scale), which map to 0 and 1
nbins = 100
count, bins, ignored = plt.hist(bkg, bins=nbins, histtype='stepfilled', edgecolor='none')
如果现在给你这些箱子的中心和计数，
xdata = .5*(bins[1:]+bins[:-1])
ydata = count
并要求您找到适合此数据的功率分布函数的参数（-> 有人告诉过您，您相信该来源），然后您可以按以下方式进行。
首先，观察功率分布函数P(x,a) 是一个单调递增的函数（即P(x1, a ) < P(x2, a) 时0 <= x1 < x2 <= 1）。这意味着上面给出的数据集已经从左到右翻转了，或者它用factor < 0 表示factor*P(x, a )。
接下来，请注意给定的数据不是在区间 [0,1] 上给出的，这对于概率密度函数来说是典型的。这意味着在尝试将幂函数分布拟合到它之前，您应该将给定的 xdata 重新缩放到 [0,1] 区间。只需通过观察图表，您就会发现 0 和 1 映射到的值是 100 和 240。但是，这只是运气，因为 matplotlib 选择了一个合理的绘图范围。当您遇到实际上不知道 0 和 1 映射到的限制时，您可以选择 xdata[0] - binwidth/2 和 xdata[-1] + binwidth/2 或（稍差的选择）xdata[0] 的不太理想（但仍然非常好的）选择和xdata[-1]。从上一段中，您知道 1 映射到 xdata[0] - binwidth/2 :=: a 和 0 映射到 xdata[-1] + binwidth/2 :=: b。执行此操作的线性映射是lambda x: (a - b)*x + b（简单代数）。
如果你将它传递给 xdata 的 [0,1] 映射版本到 curve_fit，它会给你一个很好的指数猜测。
def get_model(nobservations, binwidth, scale, offset):
def model(bin_centers, exponent):
x = (bin_centers - offset)/scale
y = exponent*x**(exponent - 1)
normed_y = nobservations * binwidth * y / np.abs(scale)
return normed_y
return model
binwidth = np.diff(xdata)[0]
p0, _ = curve_fit(get_model(nsamples, binwidth, scale=-xdata.ptp() - binwidth, offset=xdata[-1] + binwidth/2), xdata, ydata)
print(p0) # prints e.g.: 3.37117679
plt.plot(xdata, get_model(nsamples, binwidth, scale=-xdata.ptp() - binwidth, offset=xdata[-1] + binwidth/2)(xdata, *p0))
此时，您已经找到了一个相当准确的分布描述
用于生成bkg的观察结果：
f(x) = offset + scale*(exponent * x**(exponent - 1))
= (xdata[-1] + binwidth/2) + (-xdata.ptp() - binwidth)*(p0[0] * x**(p0[0] - 1))
~ 234.85 - 1.34.85*(3.37 * x**(3.37 - 1))
顺便说一下，我想指出复制bkg（来自分布的观察）
作为一个完美的副本，只有在您知道分布的确切参数（240、-140 和 3.5）并且将随机数生成的种子设置为等于初始调用之前有效的种子时，您才能做到这一点np.random.power.
如果您想使用splines 将曲线拟合到直方图，您应该从生成的样条中检索节点和系数，并将它们传递给bspleval 的函数，如here 所示。然而，写出这些方程式的主题很长，互联网上有许多资源可供您查看以了解它是如何完成的。不用说，bspleval 功能是您需要的，以防您想走那条路。如果是我，我会走上图的曲线拟合路线。
【讨论】：