如何在Python或R中找到半范数分布的第三个标准差？

Question

我有一个数据框：

import pandas as pd
import matplotlib.pyplot as plt 

df = pd.DataFrame([100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, 1, 0.9, 0.07, 0.001, 0.00001], columns=['Frequency'])

我假设该数据遵循半范数分布。

df.sort_values(by='Frequency').plot.bar(y='Frequency')
plt.savefig('test_dist.png')

此可视化结果：

如何首先对该数据拟合半范数（半钟形）分布，然后返回第三个标准差？我可以找到

scipy.stats.halfnorm

，但不知道如何使用它。如果您知道 R 的解决方案，我愿意接受其他语言，例如 R。

Answer 1

如何首先在此数据上拟合半范数（半钟形）分布

SciPy 具有将分布拟合到原始样本的函数；它不接受直方图或加权样本。假设观测值是整数 0-15（如图所示）并且某种近似值是可以的，我们可以将每个观测值乘以最接近的整数频率。

import numpy as np
freqs = [100, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, 1, 0.9, 0.07, 0.001, 0.00001]
freqs = np.asarray(freqs).round().astype(int)  # approximation: round to integers
# freqs = (np.asarray(freqs)*1000).astype(int)  # less approximation, slow
# freqs = (np.asarray(freqs)*100000).astype(int) # no approximation, very slow
vals = np.arange(len(freqs))
sample = np.repeat(vals, freqs)

现在

sample

是一个包含 100 个 0、90 个 1、...、0 个 14 和 0 个 15 的数组。

拟合半范数分布：

from scipy import stats
# assume `loc` is fixed at 0
loc, scale = stats.halfnorm(sample, floc=0)

我将把“第三个标准差”解释为相应（完全）正态分布的标准差的三倍。

scale

参数是相应正态分布的标准差，因此您需要该值的三倍：

3 * scale  # 12.069124831118353

为了确保我们在同一页面上，如下所示：

import matplotlib.pyplot as plt
dist = stats.halfnorm(loc, scale)
pdf = dist.pdf(vals)
plt.plot(vals, pdf, label='fitted half-normal')
plt.hist(sample, bins=vals, density=True, label='histogram')
plt.xlabel('value')
plt.ylabel('frequency density')
plt.title('Histogram of Data and fitted Half-Normal PDF')
plt.legend()

3 * scale

处的生存函数（右尾留下的概率质量）为：

dist.sf(3 * scale)  # 0.0026997960632601866

如果这是一个完全正态分布，那么它是 3σ 时的两倍。