python：瑞利拟合（直方图）

Question

我还在 python 学习编程。我第一次尝试使用直方图和拟合！

特别是，我有一个数据集，并制作了它的直方图。在这一点上，我应该做一个瑞利拟合，但我想不出正确设置参数的正确方法。我看loc和scale，应该是fit的参数，一般都设置为0和1，显然这样，fit效果不好！！！有人可以帮助我吗？明确地说，我附上了我正在使用的代码。

谢谢。

import os 
import numpy as np
import nrrd
import nibabel as nib
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import SimpleITK as sitk
import scipy.stats
from scipy.stats import rayleigh
import math


#fit

# Sample from this Random variable
x0 = np.array(fondi)

# 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())

# Display Results:
fig, axe = plt.subplots()

axe.hist(x0,density=1, label="background")

axe.plot(xl,stats.rayleigh(scale=scale, loc=loc).pdf(xl), label="Rayleigh")
axe.set_title("Distribution Fit")
axe.set_xlabel("Intensità")

axe.legend()
axe.grid()

我的数据 (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 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]

输出：loc=6.783540954380711 scale=6.430045149216335

Answer 1

调整MCVE

下面是一个从 Rayleigh distribution and then find its parameters using Maximum Likelihood Estimation provided by the 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 对于直方图 bin 来说是一个巨大的数字，它会降低表示的质量，因为您将拥有空的或填充较少的 bin。它还可以使统计测试（例如卡方拟合优度）由于代表性不足而失败。你当然可以让 matplotlib 估计 bins 的数量；
分布通常采用位置和比例参数，在 scipy.stats 中，他们尽最大努力以这种方式规范化 - 如果可能的话 - 每个可用的分布。要找出与usual parametric distribution definition的对应关系，需要解决以下问题：pdf(x) = pdf(y)/scale where y = (x-loc)/scale。在这种情况下，您会看到 scale 参数等同于 sigma 并且这对于原点偏移是不变的（不依赖于 loc 值）；
要调整分布，您需要在某个点执行一些 analytical/statistical 过程以从采样数据中估计参数。您的代码中缺少此部分（请参阅上面 MCVE 中的 stats.rayleigh.fit(x0)）。这部分独立于 matplotlib 绘制的任何图形，它由 scipy 处理，后者对完整数据集执行 MLE（这就是为什么更改 bin 只会影响直方图显示而不会影响其他）。

更新

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

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)

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

让我们通过统计测试来确认一下。首先我们可以检查 ECDF is compatible with the CDF of the adjusted distribution using the Kolmogorov-Smirnov Test:

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

我们还可以评估调整后分布的预期计数，并使用 Chi Square Test:

将它们与 osberved 进行比较

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)

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

两个测试都非常低且相似 p-value，这意味着不太可能满足原假设（因此它告诉我们应该拒绝它们）：

Kolmogorov: H0 = 样本取自参考分布；
卡方： H0 = 类的观察分布和预期分布没有差异；

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

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

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

在那种情况下我们无法拒绝 H0，我们有信心采样数据可能来自调整后的瑞利分布。

python：瑞利拟合（直方图）

python: rayleigh fit (histogram)

python

numpy

histogram

curve-fitting

scipy

调整MCVE

备注

更新