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

不确定度分析：我花了一年时间研究不确定性估算xs ys ts x_scale t_scale = generate_time_series() def model(xs k m): return k * xs m def l2_loss(tup xs ys): k m = tup delta = model(xs k m) - ys return numpy.dot(delta delta) k_hat m_hat = scipy.optimize.minimize(l2_loss (0 0) args=(xs ys)).x pyplot.scatter(ts ys alpha=0.5 s=100) pyplot.plot(t_scale model(x_scale k_hat m_hat) color='red' linewidth=5 alpha=0.5) p

神奇吧，这些图表与之前的图表非常相似！（正如我们本该期待的那样）

Bootstrapping算法很不错，因为它可以让你回避了任何关于生成数据的概率分布的问题。虽然它可能有点慢，但它基本上是即插即用的，并且几乎适用于所有问题。

不过请注意bootstrapping也存在不可用的情形。我的理解是，当样本数量趋于无穷大时，bootstrapping会收敛到正确的估计值。但如果你使用的是小样本集，你会得到不靠谱的结果。我通常不会相信对任何少于50个样本的问题的bootstrapping再抽样，所以你最好也不要这样做。

边注，Seaborn的 barplot实际上使用bootstrapping来绘制置信区间：

seaborn.barplot(data=d x='Month' y='Weight (kg)')

同样，Seaborn非常适合进行探索性分析，其中一些图表可以用于基本的统计。

回归

让我们提高一个档次。我们将在这一群点上做一个线形回归。

我将以我认为最普遍的方式来做这件事。我们将定义一个模型（在这种情况下是一条直线），一个损失函数（与该直线的平方偏差），然后使用通用求解器（scipy.optimize.minimize）对其进行优化。

xs ys ts x_scale t_scale = generate_time_series() def model(xs k m): return k * xs m def l2_loss(tup xs ys): k m = tup delta = model(xs k m) - ys return numpy.dot(delta delta) k_hat m_hat = scipy.optimize.minimize(l2_loss (0 0) args=(xs ys)).x pyplot.scatter(ts ys alpha=0.5 s=100) pyplot.plot(t_scale model(x_scale k_hat m_hat) color='red' linewidth=5 alpha=0.5) pyplot.ylabel('Weight of elephant (kg)')

具有不确定性的线性回归，使用最大似然方法

我们只拟合k和m，但这里没有不确定性估计。有几件事我们可以估计不确定性，但让我们从预测值的不确定性开始。

我们可以通过在拟合k和m的同时在直线周围拟合正态分布来做到这一点。我将使用最大似然方法来做到这一点。如果你不熟悉这种方法，不要害怕！如果统计学存在方法能够很容易实现（它是基本概率理论）并且很有用，那就是这种方法。

实际上，最小化平方损失（我们刚刚在前面的片段中做过）实际上是最大可能性的特殊情况！最小化平方损失与最大化所有数据概率的对数是一回事。这通常称为“对数似然”。

所以我们已经有一个表达式来减少平方损失。如果我们使方差为未知变量σ2，我们可以同时拟合它！我们现在要尽量减少的数量变成了

在这里 ^yi=kxi m是我们模型的预测值。让我们尝试拟合它！

import scipy.optimize def neg_log_likelihood(tup xs ys): # Since sigma > 0 we use use log(sigma) as the parameter instead. # That way we have an unconstrained problem. k m log_sigma = tup sigma = numpy.exp(log_sigma) delta = model(xs k m) - ys return len(xs)/2*numpy.log(2*numpy.pi*sigma**2) \ numpy.dot(delta delta) / (2*sigma**2) k_hat m_hat log_sigma_hat = scipy.optimize.minimize( neg_log_likelihood (0 0 0) args=(xs ys) ).x sigma_hat = numpy.exp(log_sigma_hat) pyplot.scatter(ts ys alpha=0.5 s=100) pyplot.plot(t_scale model(x_scale k_hat m_hat) color='green' linewidth=5) pyplot.fill_between( t_scale model(x_scale k_hat m_hat) - 1.96*sigma_hat model(x_scale k_hat m_hat) 1.96*sigma_hat color='red' alpha=0.3) pyplot.legend() pyplot.ylabel('Weight of elephant (kg)')

这里的不确定性估计实际上不是100%，因为它没有考虑k，m和σ本身的不确定性。这是一个不错的近似，但要做到正确，我们需要同时做这些事情。所以我们这样做吧。

再次启用Bootstrapping！

所以让我们把它提升到一个新的水平，并尝试估计k和m和σ的不确定性估计！ 我认为这将显示bootstrapping基本上是如何切割 - 你可以将它插入几乎任何东西，以估计不确定性。

对于每个bootstrap估计，我将绘制一条线。我们也可以采用所有这些线并计算置信区间：

pyplot.scatter(ts ys alpha=0.5 s=100) xys = list(zip(xs ys)) curves = [] for i in range(100): # sample with replacement bootstrap = [random.choice(xys) for _ in xys] xs_bootstrap = numpy.array([x for x y in bootstrap]) ys_bootstrap = numpy.array([y for x y in bootstrap]) k_hat m_hat = scipy.optimize.minimize( l2_loss (0 0) args=(xs_bootstrap ys_bootstrap) ).x curves.append(model(x_scale k_hat m_hat)) # Plot individual lines for curve in curves: pyplot.plot(t_scale curve alpha=0.1 linewidth=3 color='green') # Plot 95% confidence interval lo hi = numpy.percentile(curves (2.5 97.5) axis=0) pyplot.fill_between(t_scale lo hi color='red' alpha=0.5) pyplot.ylabel('Weight of elephant (kg)')