4

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

在下面的代码中,我生成了 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')

这是带有红色拟合函数的直方图。

之后,我想使用卡方检验来测试这个函数与直方图的拟合程度。此测试使用这些点中的观察值和预期值。为了计算期望值,我首先计算每个 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 的卡方函数中,我得到的 p 值大约为 e-5 到 e-15 数量级,这告诉我拟合函数不描述直方图:

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

但这不是真的,函数非常适合直方图!

有人知道我在哪里做错了吗?

非常感谢!!查尔斯

ps:我把delta自由度的个数设置为2,因为2个参数a1和b1是从样本中估计出来的。我尝试使用其他 ddof,但结果仍然很差!

4

1 回答 1

7

您对数组端点的计算减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)
于 2017-03-19T20:07:59.727 回答