不确定度分析：我花了一年时间研究不确定性估算

不确定度分析：我花了一年时间研究不确定性估算让我们看看我经常使用的一种数据集：转化。为了论证，我们假设正在进行一个有一定影响的A / B测试，并且我们正试图了解各州对转化率的影响。转化结果是0或1。生成此数据集的代码并不是非常重要，因此不要过多关注：所有结果为0或1时的置信区间顺便说一句 – 数值1.96是怎么来的？它与不确定性估计的大小直接相关。± 1.96意味着你将覆盖概率分布的95％左右。def plot_confidence_interval(observations_by_group): groups = list(sorted(observations_by_group.keys())) lo_bound = [] hi_bound = [] for group in groups: series = observations_by_group[group] mu std n = numpy.mean(s

最后的图表显示了数据集的分布。现在让我们试着弄清楚一个非常常见的估算器的不确定性：均值！

计算均值的不确定性 - 正态分布

在一些宽松的假设下（我一会儿回来仔细研究它），我们可以计算均值估计量的置信区间：

这里¯ X是平均值和σ是标准差，也就是方差的平方根。我不认为记住这个公式非常重要，但我觉得记住置信区间的大小与样本数的平方根成反比这个关系还是有点用的。例如，这在当你运行的A/B测试是有用的-如果你要检测1％的差异，那么你需要大约0.01−2=10 000个样本。（这只是经验值，不要把它用于你的医疗设备软件）。

顺便说一句 – 数值1.96是怎么来的？它与不确定性估计的大小直接相关。± 1.96意味着你将覆盖概率分布的95％左右。

def plot_confidence_interval(observations_by_group): groups = list(sorted(observations_by_group.keys())) lo_bound = [] hi_bound = [] for group in groups: series = observations_by_group[group] mu std n = numpy.mean(series) numpy.std(series) len(series) lo_bound.append(mu - 1.96*std*n**-0.5) hi_bound.append(mu 1.96*std*n**-0.5) pyplot.fill_between(groups lo_bound hi_bound alpha=0.2 label='Confidence interval (normal)') pyplot.scatter(ts ys alpha=0.5 s=100) observations_by_month = {} for month y in zip(d['Month'] d['Weight (kg)']): observations_by_month.setdefault(month []).append(y) plot_confidence_interval(observations_by_month) pyplot.ylabel('Weight of elephant (kg)') pyplot.legend()

请注意，这是指均值的不确定性，这与数据分布本身不是一回事。这就是为什么你看到在红色阴影区域内的蓝色点数远少于95％。如果我们添加更多的点，红色阴影区域将变得越来越窄，而其中蓝色点数仍将具有差不多的比例。然而，理论上真正的平均值应该有95％时间处于红色阴影区域内。

我之前提到，置信区间的公式仅适用于一些宽松的假设。那些假设是什么？是正态的假设。根据中心极限定理，这对于大量的观测值也是可行的。

所有结果为0或1时的置信区间

让我们看看我经常使用的一种数据集：转化。为了论证，我们假设正在进行一个有一定影响的A / B测试，并且我们正试图了解各州对转化率的影响。转化结果是0或1。生成此数据集的代码并不是非常重要，因此不要过多关注：

STATES = ['CA' 'NY' 'FL' 'TX' 'PA' 'IL' 'OH'] GROUPS = ['test' 'control'] def generate_binary_categorical(states=STATES groups=GROUPS k=400 zs=[0 0.2] z_std=0.1 b=-3 b_std=1): # Don't pay too much attention to this code. The main thing happens in # numpy.random.binomial which is where we draw the "k out of n" outcomes. output = {} e_obs_per_state = numpy.random.exponential(k size=len(states)) state_biases = numpy.random.normal(b b_std size=len(states)) for group z in zip(groups zs): noise = numpy.random.normal(z z_std size=len(states)) ps = 1 / (1 numpy.exp(-(state_biases noise))) ns = numpy.random.poisson(e_obs_per_state) ks = numpy.random.binomial(ns ps) output[group] = (ns ks) return output

对于每个州和每个“组”（测试和控制），我们生成了n个用户，其中k个已转化。让我们绘制每个州的转化率，看看发生了什么！

data = generate_binary_categorical() for group (ns ks) in data.items(): pyplot.scatter(STATES ks/ns label=group alpha=0.7 s=400) pyplot.ylabel('Conversion rate') pyplot.legend()

由于所有结果都是0或1，并且以相同（未知）概率绘制，我们知道1和0的数量遵循二项分布。这意味着“n个用户中 k个已转化”的情形的置信区间是Beta分布。

记住置信区间的公式使我获益良多，而且我觉得比起我以前用的（基于法线的）公式，我可能更倾向用它。特别需要记住的是：

n k = 100 3 scipy.stats.beta.ppf([0.025 0.975] k n-k) array([0.00629335 0.07107612])

代入n和k的值可以算出95％的置信区间。在这个例子里，我们看到如果我们有100个网站访问者，其中3个购买了产品，那么置信区间就是0.6％-7.1％。

让我们用我们的数据集试试：

for group (ns ks) in data.items(): lo = scipy.stats.beta.ppf(0.025 ks ns-ks) hi = scipy.stats.beta.ppf(0.975 ks ns-ks) mean = ks/ns pyplot.errorbar(STATES y=mean yerr=[mean-lo hi-mean] label=group alpha=0.7 linewidth=0 elinewidth=50) pyplot.ylabel('Conversion rate') pyplot.legend()

哎哟不错噢~

Bootstrapping算法（拔靴法）

另一种有用的方法是bootstrapping算法（拔靴法），它允许你在无需记忆任何公式的情况下做同样的统计。这个算法的核心是计算均值，但是是为n次再抽样（bootstrap）计算均值，其中每个bootstrap是我们观测中的随机样本（替换）。对于每次自助抽样，我们计算平均值，然后将在97.5%和2.5%之间的均值作为置信区间：

lo_bound = [] hi_bound = [] months = sorted(observations_by_month.keys()) for month in months: series = observations_by_month[month] bootstrapped_means = [] for i in range(1000): # sample with replacement bootstrap = [random.choice(series) for _ in series] bootstrapped_means.append(numpy.mean(bootstrap)) lo_bound.append(numpy.percentile(bootstrapped_means 2.5)) hi_bound.append(numpy.percentile(bootstrapped_means 97.5)) pyplot.scatter(ts ys alpha=0.5 s=100) pyplot.fill_between(months lo_bound hi_bound alpha=0.2 label='Confidence interval') pyplot.ylabel('Weight of elephant (kg)') pyplot.legend()