Python 中的卡方拟合优度检验:p 值太低,但拟合函数正确
Chi-squared goodness of fit test in Python: way too low p-values, but the fitting function is correct
尽管在相关问题中搜索了两天,但我还没有真正找到这个问题的答案...
在下面的代码中,我生成了 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的正态分布,并以scaling_factor进行缩放以满足样本未归一化的事实。
它确实非常符合直方图:
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')
Here's the plot of the histogram with the fitting function in red.
之后,我想使用卡方检验来测试此函数与直方图的拟合程度。
此测试使用这些点中的观察值和预期值。为了计算期望值,我首先计算每个 bin 的中间位置,这个信息包含在数组 x_middle 中。然后我计算每个 bin 中间点的拟合函数值,它给出 expected_value 数组:
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))
但事实并非如此,该函数非常符合直方图!
有人知道我哪里弄错了吗?
非常感谢!!
查尔斯
p.s.:我把delta自由度数设置为2,因为2个参数a1和b1是从样本中估计出来的。我尝试使用其他ddof,但结果仍然很差!
你对数组终点的计算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)
尽管在相关问题中搜索了两天,但我还没有真正找到这个问题的答案...
在下面的代码中,我生成了 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的正态分布,并以scaling_factor进行缩放以满足样本未归一化的事实。 它确实非常符合直方图:
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')
Here's the plot of the histogram with the fitting function in red.
之后,我想使用卡方检验来测试此函数与直方图的拟合程度。 此测试使用这些点中的观察值和预期值。为了计算期望值,我首先计算每个 bin 的中间位置,这个信息包含在数组 x_middle 中。然后我计算每个 bin 中间点的拟合函数值,它给出 expected_value 数组:
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))
但事实并非如此,该函数非常符合直方图!
有人知道我哪里弄错了吗?
非常感谢!! 查尔斯
p.s.:我把delta自由度数设置为2,因为2个参数a1和b1是从样本中估计出来的。我尝试使用其他ddof,但结果仍然很差!
你对数组终点的计算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)