如何使用 Astropy Kuiper 测试
How to use Astropy Kuiper test
我想测试数据集的分布与均值 = 0 和方差 = 1 的高斯分布的接近程度。
来自 astropy.stats
的 kuiper 测试有一个 cdf
参数,来自文档:“一个可调用的,用于评估正在测试的分布的 CDF。将使用一次所有值的向量。默认是均匀分布”,但我不知道如何使用它来测试正态分布。如果我想要例如均值为 0.2,方差为 2 的正态分布?
所以我使用了 kuiper_two,也是来自 astropy,并创建了一个随机正态分布。请参阅下面的示例。
我看到的问题是它取决于我生成的要比较的数据点的数量。如果我使用 100 个而不是 10000 个数据点,则概率 (fpp) 会提高到 43%。
我想问题是,我该如何正确地做到这一点?另外,我如何解释 D 号?
# create data and its cdf
np.random.seed(0)
data = np.random.normal(loc=0.02, scale=1.15, size=50)
data_sort = np.sort(data)
data_cdf = [x/len(data) for x in range(0, len(data))]
# create the normal data with mean 0 and variance 1
xx = np.random.normal(loc=0, scale=1, size=10000)
xx_sort = np.sort(xx)
xx_cdf = [x/len(xx) for x in range(0, len(xx))]
# compute the pdf for a plot
x = np.linspace(-4, 4, 50)
x_pdf = stats.norm.pdf(x, 0, 1)
# we can see it all in a plot
fig, ax = plt.subplots(figsize=(8, 6))
plt.hist(xx, bins=20, density=True, stacked=True, histtype='stepfilled', alpha=0.6)
plt.hist(data, density=True, stacked=True, histtype='step', lw=3)
plt.plot(x, x_pdf, lw=3, label='G($\mu=0$, $\sigma^2=1$)')
ax2 = ax.twinx()
ax2.plot(xx_sort, xx_cdf, marker='o', ms=8, mec='green', mfc='green', ls='None')
ax2.plot(data_sort, data_cdf, marker='^', ms=8, mec='orange', mfc='orange', ls='None')
# Kuiper test
D, fpp = kuiper_two(data_sort, xx_sort)
print('# D number =', round(D, 5))
print('# fpp =', round(fpp, 5))
# Which resulted in:
# D number = 0.211
# fpp = 0.14802
astropy.stats.kuiper
期望第一个参数是要测试的分布样本,第二个参数是要测试的分布的 CDF。
此变量是一个 Callable
并期望其自身 (a) 个样本 return 累积分布函数下的值。您可以为此使用 scipy.stat
的 CDF。通过使用 functools.partial
我们可以设置任何参数。
from scipy import stats
from scipy.stats import norm
from astropy.stats import kuiper
from functools import partial
from random import shuffle
np.random.seed(0)
data = np.random.normal(loc=0.02, scale=1.15, size=50)
print(kuiper(data, partial(norm.cdf, loc=0.2, scale=2.0)))
# Output: (0.2252118027033838, 0.08776036566607946)
# The data does not have to be sorted, in case you wondered:
shuffle(data)
print(kuiper(data, partial(norm.cdf, loc=0.2, scale=2.0)))
# Output: (0.2252118027033838, 0.08776036566607946)
在这张来自 the Wikipedia article about this test 的图表中,您可以了解 Kuiper 统计量 V 测量的内容:
[]3
如果共享比较分布的参数,则距离会更小,并且各个基础 CDF 相同的估计概率会增加:
print(kuiper(data, partial(norm.cdf, loc=0.02, scale=1.15)))
# Output: (0.14926352419821276, 0.68365004302431)
函数astropy.stats.kuiper_two
相反,需要两个经验数据样本相互比较。因此,如果您想与具有易于处理的 CDF 的分布进行比较,最好直接使用 CDF(使用 kuiper
)而不是从比较分布中采样(并使用 kuiper_two
)。
还有一个吹毛求疵。除了在不是 CDF 的变量名中使用 CDF 之外,这个公式比上面的公式更具可读性:
data_cdf = np.linspace(0.0, 1.0, len(data), endpoint=False)
xx_cdf = np.linspace(0.0, 1.0, len(xx), endpoint=False)
我想测试数据集的分布与均值 = 0 和方差 = 1 的高斯分布的接近程度。
来自 astropy.stats
的 kuiper 测试有一个 cdf
参数,来自文档:“一个可调用的,用于评估正在测试的分布的 CDF。将使用一次所有值的向量。默认是均匀分布”,但我不知道如何使用它来测试正态分布。如果我想要例如均值为 0.2,方差为 2 的正态分布?
所以我使用了 kuiper_two,也是来自 astropy,并创建了一个随机正态分布。请参阅下面的示例。
我看到的问题是它取决于我生成的要比较的数据点的数量。如果我使用 100 个而不是 10000 个数据点,则概率 (fpp) 会提高到 43%。
我想问题是,我该如何正确地做到这一点?另外,我如何解释 D 号?
# create data and its cdf
np.random.seed(0)
data = np.random.normal(loc=0.02, scale=1.15, size=50)
data_sort = np.sort(data)
data_cdf = [x/len(data) for x in range(0, len(data))]
# create the normal data with mean 0 and variance 1
xx = np.random.normal(loc=0, scale=1, size=10000)
xx_sort = np.sort(xx)
xx_cdf = [x/len(xx) for x in range(0, len(xx))]
# compute the pdf for a plot
x = np.linspace(-4, 4, 50)
x_pdf = stats.norm.pdf(x, 0, 1)
# we can see it all in a plot
fig, ax = plt.subplots(figsize=(8, 6))
plt.hist(xx, bins=20, density=True, stacked=True, histtype='stepfilled', alpha=0.6)
plt.hist(data, density=True, stacked=True, histtype='step', lw=3)
plt.plot(x, x_pdf, lw=3, label='G($\mu=0$, $\sigma^2=1$)')
ax2 = ax.twinx()
ax2.plot(xx_sort, xx_cdf, marker='o', ms=8, mec='green', mfc='green', ls='None')
ax2.plot(data_sort, data_cdf, marker='^', ms=8, mec='orange', mfc='orange', ls='None')
# Kuiper test
D, fpp = kuiper_two(data_sort, xx_sort)
print('# D number =', round(D, 5))
print('# fpp =', round(fpp, 5))
# Which resulted in:
# D number = 0.211
# fpp = 0.14802
astropy.stats.kuiper
期望第一个参数是要测试的分布样本,第二个参数是要测试的分布的 CDF。
此变量是一个 Callable
并期望其自身 (a) 个样本 return 累积分布函数下的值。您可以为此使用 scipy.stat
的 CDF。通过使用 functools.partial
我们可以设置任何参数。
from scipy import stats
from scipy.stats import norm
from astropy.stats import kuiper
from functools import partial
from random import shuffle
np.random.seed(0)
data = np.random.normal(loc=0.02, scale=1.15, size=50)
print(kuiper(data, partial(norm.cdf, loc=0.2, scale=2.0)))
# Output: (0.2252118027033838, 0.08776036566607946)
# The data does not have to be sorted, in case you wondered:
shuffle(data)
print(kuiper(data, partial(norm.cdf, loc=0.2, scale=2.0)))
# Output: (0.2252118027033838, 0.08776036566607946)
在这张来自 the Wikipedia article about this test 的图表中,您可以了解 Kuiper 统计量 V 测量的内容:
[
如果共享比较分布的参数,则距离会更小,并且各个基础 CDF 相同的估计概率会增加:
print(kuiper(data, partial(norm.cdf, loc=0.02, scale=1.15)))
# Output: (0.14926352419821276, 0.68365004302431)
函数astropy.stats.kuiper_two
相反,需要两个经验数据样本相互比较。因此,如果您想与具有易于处理的 CDF 的分布进行比较,最好直接使用 CDF(使用 kuiper
)而不是从比较分布中采样(并使用 kuiper_two
)。
还有一个吹毛求疵。除了在不是 CDF 的变量名中使用 CDF 之外,这个公式比上面的公式更具可读性:
data_cdf = np.linspace(0.0, 1.0, len(data), endpoint=False)
xx_cdf = np.linspace(0.0, 1.0, len(xx), endpoint=False)