尽管我花了两天时间研究了相关的问题,但我还没有找到这个问题的答案…
在下面的代码中,我生成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')

Here's the plot of the histogram with the fitting function in red.
之后,我想用卡方检验来检验这个函数与直方图的拟合程度。
此测试使用这些点的观测值和期望值。为了计算期望值,我首先计算每个箱子中间的位置,这个信息包含在数组x_middle中。然后,在每个bin的中间点计算拟合函数的值,这给出了预期值数组:
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)

把这个插入到函数的函数中,我得到p值大约是1到3个数量级,这告诉我拟合函数没有描述直方图:
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)

注意- 1的第二个参数中的额外linspace()
更简洁的版本是
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)

09-30 15:00
查看更多