python：瑞利拟合（直方图）

Question

我\xe2\x80\x99m 仍在尝试用Python 编程。\n我第一次尝试使用直方图和拟合！

\n

特别是，我有一个数据集，并制作了它的直方图。此时我应该进行瑞利拟合，但我无法找出正确设置参数的正确方法。我读到loc和scale，这应该是fit的参数，通常设置为0和1。显然，这样一来，fit\xe2\x80\x99就不能很好地工作了！有人可以帮助我吗？\n为了明确起见，我附上了我正在使用的代码。

\n

谢谢。

\n

import os \nimport numpy as np\nimport nrrd\nimport nibabel as nib\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nimport SimpleITK as sitk\nimport scipy.stats\nfrom scipy.stats import rayleigh\nimport math\n\n\n#fit\n\n# Sample from this Random variable\nx0 = np.array(fondi)\n\n# Adjust Distribution parameters\nloc, scale = stats.rayleigh.fit(x0) # (9.990726961181025, \n4.9743913760956335)\n\n# Tabulate over sample range (PDF display):\nxl = np.linspace(x0.min(), x0.max())\n\n# Display Results:\nfig, axe = plt.subplots()\n\naxe.hist(x0,density=1, label="background")\n\naxe.plot(xl,stats.rayleigh(scale=scale, loc=loc).pdf(xl), label="Rayleigh")\naxe.set_title("Distribution Fit")\naxe.set_xlabel("Intensit\xc3\xa0")\n\naxe.legend()\naxe.grid()\n

\n

我的数据（fondi）是这样的：[13 15 13 14 12 13 12 14 15 12 11 10 11 15 18 11 11 11 13 15 15 15 11 12\n13 12 15 15 15 12 12 11 14 16 11 13 14 1 6 17 24 21 16 20 18 18 19 21 22\n19 15 16 15 13 14 16 18 21 19 22 14 13 14 15 14 17 19 17 16 18 12 15 17\n17 16 17 16 19 17 14 13 1 6 16 13 15 17 17 20 18 17 12 19 14 15 15 14 13\n17 16 14 12 11 12 20 19 16 24 19 20 19 17 16 17 16 19 22 17 16 20 22 21\n22 20 14 18 16 19 20 17 2 0 22 20 22 19 17 13 16 18 14 16 20 20 18 19 19\n16 19 12 12 14 14 13 15 16 16 19 16 17 12 11 11 10 12 11 11 13 14 13 17\n8 8 8 10 10 10 14 16 11 9 9 11 10 17 13 15 19 15 13 16 17 14 12 13\n14 11 10 15 13 12 12 11 10 9 9 9 9 8 15 16 12 9 11 9 10 10 7 7\n7 21 19 13 10 15 12 10 10 9 8 10 2 0 14 13 11 13 15 14 10 11 12 16 17\n15 12 13 16 15 13 14 17 14 13 15 13 11 14 15 17 18 22 21 16 17 22 17 17\n18 26 17 19 21 16 15 19 1 9 22 19 18 17 18 18 12 17 17 17 18 14 16 20 17\n16 16 18 16 19 18 18 20 18]

\n

输出：loc=6.783540954380711 比例=6.430045149216335

\n

\n

Answer 1

调整MCVE

下面是从瑞利分布中绘制试验数据集，然后使用该方法提供的最大似然估计scipy.stats.rv_continuous.fit找到其参数的简单过程：

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

# Create a Continuous Variable: 
X = stats.rayleigh(loc=10, scale=5)

# Sample from this Random variable
x0 = X.rvs(size=10000, random_state=123)

# Adjust Distribution parameters
loc, scale = stats.rayleigh.fit(x0) # (9.990726961181025, 4.9743913760956335)

# Tabulate over sample range (PDF display):
xl = np.linspace(x0.min(), x0.max(), 100)

# Display Results:
fig, axe = plt.subplots()
axe.hist(x0, density=1, label="Sample")
axe.plot(xl, X.pdf(xl), label="Exact Distribution")
axe.plot(xl, stats.rayleigh(scale=scale, loc=loc).pdf(xl), label="Adjusted Distribution")
axe.set_title("Distribution Fit")
axe.set_xlabel("Variable, $x$ $[\mathrm{AU}]$")
axe.set_ylabel("Density, $f(x)$ $[\mathrm{AU}^{-1}]$")
axe.legend()
axe.grid()

它呈现如下：

笔记

我想提请大家注意以下几个要点：

• 对于直方图箱来说，300 是一个巨大的数字，它会降低表示的质量，因为你将拥有空的或填充的箱。它还可能导致统计检验（例如卡方拟合优度）因箱的代表性不足而失败。您当然可以估计matplotlib垃圾箱的数量；
• 分布通常采用位置和尺度参数，在scipy.stats可能的情况下，它们会尽力以这种方式标准化每个可用的分布。为了找出 与通常参数分布定义的对应关系，需要解决以下问题：pdf(x) = pdf(y)/scale其中y = (x-loc)/scale。在这种情况下，您将看到该scale参数等效于sigma并且这对于原点移位是不变的（不依赖于该loc值）；
• 要调整分布，您需要在某个时刻执行一些分析/统计程序，以根据采样数据估计参数。您的代码中缺少此部分（请参阅stats.rayleigh.fit(x0)上面的 MCVE）。这部分独立于 绘制的任何图形matplotlib，它是由scipy对完整数据集执行 MLE 来处理的（这就是为什么更改 bin 只影响直方图显示而不影响其他内容）。

更新

根据您的帖子更新，我完成了我的答案。使用您提供的数据集：

x0 = np.array([13, 15, 13, 14, 12, 13, 12, 14, 15, 12, 11, 10, 11, 15, 18, 11, 11, 11, 13,
               15, 15, 15, 11, 12, 13, 12, 15, 15, 15, 12, 12, 11, 14, 16, 11, 13, 14, 16,
               17, 24, 21, 16, 20, 18, 18, 19, 21, 22, 19, 15, 16, 15, 13, 14, 16, 18, 21,
               19, 22, 14, 13, 14, 15, 14, 17, 19, 17, 16, 18, 12, 15, 17, 17, 16, 17, 16,
               19, 17, 14, 13, 16, 16, 13, 15, 17, 17, 20, 18, 17, 12, 19, 14, 15, 15, 14,
               13, 17, 16, 14, 12, 11, 12, 20, 19, 16, 24, 19, 20, 19, 17, 16, 17, 16, 19,
               22, 17, 16, 20, 22, 21, 22, 20, 14, 18, 16, 19, 20, 17, 20, 22, 20, 22, 19,
               17, 13, 16, 18, 14, 16, 20, 20, 18, 19, 19, 16, 19, 12, 12, 14, 14, 13, 15,
               16, 16, 19, 16, 17, 12, 11, 11, 10, 12, 11, 11, 13, 14, 13, 17, 8, 8, 8, 10,
               10, 10, 14, 16, 11, 9, 9, 11, 10, 17, 13, 15, 19, 15, 13, 16, 17, 14, 12, 13,
               14, 11, 10, 15, 13, 12, 12, 11, 10, 9, 9, 9, 9, 8, 15, 16, 12, 9, 11, 9, 10,
               10, 7, 7, 7, 21, 19, 13, 10, 15, 12, 10, 10, 9, 8, 10, 20, 14, 13, 11, 13, 15,
               14, 10, 11, 12, 16, 17, 15, 12, 13, 16, 15, 13, 14, 17, 14, 13, 15, 13, 11, 14,
               15, 17, 18, 22, 21, 16, 17, 22, 17, 17, 18, 26, 17, 19, 21, 16, 15, 19, 19, 22,
               19, 18, 17, 18, 18, 12, 17, 17, 17, 18, 14, 16, 20, 17, 16, 16, 18, 16, 19, 18,
               18, 20, 18])

我们可以尝试调整瑞利分布：

p = stats.rayleigh.fit(x0)
X = stats.rayleigh(*p)

从视觉上看，贴合度不太好：

让我们通过统计测试来证实这一点。首先，我们可以使用Kolmogorov-Smirnov测试检查ECDF是否与调整后的分布的 CDF 兼容：

kst = stats.kstest(x0, X.cdf)
# KstestResult(statistic=0.12701044409231593, pvalue=0.0001232197856051324)

我们还可以评估调整后分布的预期计数，并使用卡方 检验将它们与观察到的计数进行比较：

c, b = np.histogram(x0)
ct = np.diff(X.cdf(b))*np.sum(c)
c2t = stats.chisquare(c, ct, ddof=2)
# Power_divergenceResult(statistic=31.874916914227434, pvalue=4.284273564311872e-05)

自由度差等于 2，因为除了卡方统计之外，我们还必须估计瑞利分布的loc和参数（因此在测试调用中）。scaleddof=2

两个检验都有非常低且相似的p 值，这意味着满足零假设的可能性很小（因此它告诉我们应该拒绝它们）：

• 柯尔莫哥洛夫： H0 = 样本是从参考分布中抽取的；
• 卡方： H0 = 观察到的分布和预期分布中的类别之间没有差异；

很难相信您的数据集来自调整后的瑞利分布。

您可以将这些结果与 MCVE 中绘制的合成数据进行比较，测试返回高于 10% 的 p 值：

# KstestResult(statistic=0.0097140857969642, pvalue=0.3019167138216704)
# Power_divergenceResult(statistic=11.170065854104491, pvalue=0.13137094282775724)

如果这种情况我们不能拒绝 H0，我们相信采样数据可能来自调整后的瑞利分布。