Python中的卡方拟合优度检验：p值太低，但拟合函数是正确的。

Question

虽然在相关问题中搜索了两天，但我还没有真正找到这个问题的答案...

在下面的代码中，我生成了 n 个正态分布的随机变量，然后将它们表示成直方图：

import numpy as np
import matplotlib.pyplot as plt

n = 10000                        # number of generated random variables 
x = np.random.normal(0,1,n)      # generate n random variables

# plot this in a non-normalized histogram:
plt.hist(x, bins='auto', normed=False)    

# get the arrays containing the bin counts and the bin edges:
histo, bin_edges = np.histogram(x, bins='auto', normed=False)
number_of_bins = len(bin_edges)-1

之后，找到了一条曲线拟合函数及其参数。该函数服从正态分布，具有参数a1和b1，并通过缩放因子进行缩放以满足未正则化的样本事实。它确实很好地拟合了直方图：

import scipy as sp

a1, b1 = sp.stats.norm.fit(x)

scaling_factor = n*(x.max()-x.min())/number_of_bins

plt.plot(x_achse,scaling_factor*sp.stats.norm.pdf(x_achse,a1,b1),'b')

这是带有红色拟合函数的直方图绘图。 之后，我想测试这个函数对直方图的拟合效果如何，使用卡方检验。该检验使用观察值和这些点上的期望值。为了计算期望值，我首先计算每个箱子中间位置的位置，这些信息包含在x_middle数组中。然后，我计算每个箱子中间点处拟合函数的值，得到期望值数组。

observed_values = histo

bin_width = bin_edges[1] - bin_edges[0]

# array containing the middle point of each bin:
x_middle = np.linspace(  bin_edges[0] + 0.5*bin_width,    
           bin_edges[0] + (0.5 + number_of_bins)*bin_width,
           num = number_of_bins) 

expected_values = scaling_factor*sp.stats.norm.pdf(x_middle,a1,b1)

将此插入Scipy的卡方函数中，我得到大约为e-5至e-15数量级的p值，这告诉我拟合函数不能描述直方图。

print(sp.stats.chisquare(observed_values,expected_values,ddof=2)) 

但事实并非如此，该函数非常适合直方图！

有人知道我错在哪里了吗？

非常感谢！ Charles

p.s.：我将自由度增量的数量设置为2，因为参数a1和b1是从样本中估计出来的。 我尝试使用其他ddof，但结果仍然很差！

- Charles M.

Answer 1

你对数组 x_middle 的终点计算有误，正确应该是：

x_middle = np.linspace(bin_edges[0] + 0.5*bin_width,    
                       bin_edges[0] + (0.5 + number_of_bins - 1)*bin_width,
                       num=number_of_bins)

注意在 linspace() 的第二个参数中有额外的 -1。

更简洁的版本是：

x_middle = 0.5*(bin_edges[1:] + bin_edges[:-1])

计算expected_values的另一种（可能更准确的）方法是使用CDF的差异，而不是使用每个间隔中间的PDF来近似这些差异：

In [75]: from scipy import stats

In [76]: cdf = stats.norm.cdf(bin_edges, a1, b1)

In [77]: expected_values = n * np.diff(cdf)

通过这个计算，我可以得到卡方检验的以下结果：

In [85]: stats.chisquare(observed_values, expected_values, ddof=2)
Out[85]: Power_divergenceResult(statistic=61.168393496775181, pvalue=0.36292223875686402)

- Warren Weckesser

---

你对于计算“expected_values”的建议确实更好。从数学角度来看，这个新得到的数组对应于拟合曲线下的积分（如果乘以n），而直方图计数应该对应于它。谢谢！ - Charles M.