不确定度分析：我花了一年时间研究不确定性估算
不确定度分析：我花了一年时间研究不确定性估算我将切换到一些贝叶斯方法，我们通过绘制样本来估计k，m和σ。它类似于bootstrapping，但MCMC有更好的理论基础（我们使用贝叶斯规则从“后验分布”中抽样），并且它通常要快几个数量级。现在它会变得有点狂野了~~这里的技巧是，对于（k，m，σ）的每个bootstrap估计，我们还需要绘制随机预测。正如您在代码中看到的，我们实际上是将随机正态变量添加到y的预测值。这也是为什么这个形状最终变成一个大波浪形的原因。不幸的是，bootstrapping对于这个问题来说相当缓慢 - 对于每个bootstrap，我们需要拟合一个模型。让我们看看另一个选项：马尔可夫链蒙特卡罗方法
哇，这里发生了什么？这种不确定性与之前的情节非常不同。这看起来很混乱，直到你意识到他们展示了两个非常不同的东西：
- 第一个图找到k和m的一个解，并显示预测的不确定性。所以，如果你被问到下个月大象体重的范围是什么，你可以从图表中得到它。
- 第二个图找到了k和m的许多解，并显示了kx m的不确定性。因此，这回答了一个不同的问题 - 大象体重随时间变化的趋势是什么，趋势的不确定性是什么？
事实证明，我们可以将这两种方法结合起来，并通过拟合绘制bootstrapping样本并同时拟合k，m和σ使其更加复杂。然后，对于每个估算，我们可以预测新值y。我们这样做：
pyplot.scatter(ts ys alpha=0.5 s=100) xys = list(zip(xs ys)) curves = [] for i in range(4000): # 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 log_sigma_hat = scipy.optimize.minimize( neg_log_likelihood (0 0 0) args=(xs_bootstrap ys_bootstrap) ).x curves.append( model(x_scale k_hat m_hat) # Note what's going on here: we're _adding_ the random term # to the predictions! numpy.exp(log_sigma_hat) * numpy.random.normal(size=x_scale.shape) ) # 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)')
哎呦，又不错！它现在变得很像回事了 - 如果仔细观察，你会看到一个双曲线形状！
这里的技巧是，对于（k，m，σ）的每个bootstrap估计，我们还需要绘制随机预测。正如您在代码中看到的，我们实际上是将随机正态变量添加到y的预测值。这也是为什么这个形状最终变成一个大波浪形的原因。
不幸的是，bootstrapping对于这个问题来说相当缓慢 - 对于每个bootstrap，我们需要拟合一个模型。让我们看看另一个选项：
马尔可夫链蒙特卡罗方法
现在它会变得有点狂野了~~
我将切换到一些贝叶斯方法，我们通过绘制样本来估计k，m和σ。它类似于bootstrapping，但MCMC有更好的理论基础（我们使用贝叶斯规则从“后验分布”中抽样），并且它通常要快几个数量级。
为此，我们将使用一个名为emcee的库，我发现它很好用。 它所需要的只是一个Log—likelihood函数，我们之前已经定义过了！ 我们只需要用它的负数。
import emcee xs ys ts x_scale t_scale = generate_time_series() def log_likelihood(tup xs ys): return -neg_log_likelihood(tup xs ys) ndim nwalkers = 3 10 p0 = [numpy.random.rand(ndim) for i in range(nwalkers)] sampler = emcee.EnsembleSampler(nwalkers ndim log_likelihood args=[xs ys]) sampler.run_mcmc(p0 10000)
让我们绘制k和m的采样值！
# Grab the last 10 from each walker samples = sampler.chain[: -10: :].reshape((-1 ndim)) pyplot.scatter(ts ys alpha=0.5 s=100) for k m log_sigma in samples: pyplot.plot(t_scale model(x_scale k m) alpha=0.1 linewidth=3 color='green') pyplot.ylabel('Weigh of elephant (kg)')
这些方法还有一些内容 - 采样有点挑剔，需要一些人工参与才能正常工作。我不想在这里深入解释所有细节，而且我自己就是一个门外汉。但它通常比booststrapping快几个数量级，并且它还可以更好地处理数据更少的情况。
我们最终得到了k，m，σ后验分布的样本。我们可以看看这些未知数的概率分布：
# Grab slightly more samples this time samples = sampler.chain[: -500: :].reshape((-1 ndim)) k_samples m_samples log_sigma_samples = samples.T seaborn.distplot(k_samples label='k') seaborn.distplot(m_samples label='m') seaborn.distplot(numpy.exp(log_sigma_samples) label='sigma') pyplot.legend()
你可以看到围绕k = 200，m = 1000，σ= 100时的分布中心，我们原本也是如此构建它们的！这看起来挺令人放心的~
最后，我们可以使用与boostraps相同的方法绘制预测的完全不确定性：
pyplot.scatter(ts ys alpha=0.5 s=100) samples = sampler.chain[: -4000: :].reshape((-1 ndim)) curves = [] for k m log_sigma in samples: curves.append( model(x_scale k m) numpy.exp(log_sigma) * numpy.random.normal(size=x_scale.shape) ) # 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)')
这些贝叶斯方法并没有在这里结束。特别是有几个库可以解决这些问题。事实证明，如果您以更结构化的方式表达问题（而不仅仅是负对数似然函数），您可以将采样比例调整为大问题（如数千个未知参数）。对于Python来说，有PyMC3和PyStan，以及稍微更小众一点的Edward和Pyro。
结语
似乎我把你们带进了一个深坑 - 但这些方法其实可以走得更远。事实上，强迫自己估计我做的任何事情的不确定性是一个很好的体验，可以逼着自己学习很多的统计知识。
根据数据做出决策很难！但如果我们对量化不确定性更加严格，我们可能会做出更好的决策。现在要做到这一点并不容易，但我真的希望我们能够使用更便捷的工具，从而看到这些方法的普及。
本文有关的所有代码都可以在Github上可以找到：
https://github.com/erikbern/uncertainty。
相关报道：
https://erikbern.com/2018/10/08/the-hackers-guide-to-uncertainty-estimates.html
