讓我從年初的一個flag說起：
我的新年目标：我在2018年期間繪制的每一幅圖表都要包含不确定性估算
為什麼立下這個flag？因為我在各種大會上聽膩了人們争論每個月微件（widget）的數量是上升還是下降，或者微件方法X是否比微件方法Y更有效率。
而且我發現，對于幾乎任何圖表，量化不确定性都很有用，是以我也開始嘗試。
然而我很快發現，我給自己挖了個深坑。幾個月後：
現在已經是今年的第4個月，我要告訴你，估算不确定性的水還挺深。
我從未學過統計學，也沒有通過機器學習來逆向了解過它。是以我算是半路出家，在慢慢自學統計知識。今年早些時候，我還隻了解一些關于Bootstrapping算法（拔靴法）和置信區間的基本知識，但随着時間的推移，我學會了蒙特卡羅方法和逆Hessians矩陣（黑塞矩陣）等全套把戲。
這些方法很有用，我也想把這一年的經營教訓分享給大家。
從資料開始
我相信沒有具體例子是無法真正學到東西的，是以讓我們先制造一些資料。我們将生成一個假的時間系列，其日期範圍從2017-07-01至2018-07-31，比如說這個序列是一頭大象重量的觀測值。
def generate_time_series(k=200, m=1000, sigma=100, n=50,
start_date=datetime.date(2017, 7, 1)):
xs = numpy.linspace(0, 1, n, endpoint=False)
ys = [k*x + m + random.gauss(0, sigma) for x in xs]
ts = [start_date + datetime.timedelta(x)*365 for x in xs]
x_scale = numpy.linspace(-1, 2, 500) # for plotting
t_scale = [start_date + datetime.timedelta(x)*365 for x in x_scale]
return xs, ys, ts, x_scale, t_scale
xs, ys, ts, x_scale, t_scale = generate_time_series()
在開始之前，我們需要做圖來看看發生了什麼！
pyplot.scatter(ts, ys, alpha=0.5, s=100)
pyplot.xlabel('Date')
pyplot.ylabel('Weight of elephant (kg)')
首先，我們不用任何花哨的模型，我們隻是将其分解為幾個區間（bucket）并計算每個區間的平均值。不過，讓我們先停下來談談不确定性。
資料分布與不确定性
之前我一直搞不清“不确定性”的意思，但我認為搞清楚這一點非常重要。我們可以為多種不同的資料估算分布：
1. 資料本身。給定一定的時間範圍（t ，t '），在這個時間間隔内大象體重的分布是什麼？
2.某些參數的不确定性。如參數k線上性關系y = k t + m裡，或者某些估算器的不确定性，就像許多觀測值的平均值一樣。
3.預測數值的不确定性。是以，如果我們預測日期為t（可能在未來）時大象的重量是y公斤，我們想知道數量y的不确定性。
讓我們從最基本的模型開始 - 隻需在區間中分解問題。如果我們隻是想學習一些關于分布和不确定性估計的基本概念，那麼我推薦Seaborn軟體包。Seaborn通常在資料幀上運作，是以我們需要進行轉換：
d = pandas.DataFrame({'x': xs, 't': ts, 'Weight (kg)': ys})
d['Month'] = d['t'].apply(lambda t: t.strftime('%Y-%m'))
seaborn.boxplot(data=d, x='Month', y='Weight (kg)')
最後的圖表顯示了資料集的分布。現在讓我們試着弄清楚一個非常常見的估算器的不确定性：均值！
計算均值的不确定性 - 正态分布
在一些寬松的假設下（我一會兒回來仔細研究它），我們可以計算均值估計量的置信區間：
這裡¯ 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)')
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.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')
由于所有結果都是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％。
讓我們用我們的資料集試試：
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)
哎喲不錯噢~
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.fill_between(months, lo_bound, hi_bound, alpha=0.2,
label='Confidence interval')
神奇吧，這些圖表與之前的圖表非常相似！（正如我們本該期待的那樣）
Bootstrapping算法很不錯，因為它可以讓你回避了任何關于生成資料的機率分布的問題。雖然它可能有點慢，但它基本上是即插即用的，并且幾乎适用于所有問題。
不過請注意bootstrapping也存在不可用的情形。我的了解是，當樣本數量趨于無窮大時，bootstrapping會收斂到正确的估計值。但如果你使用的是小樣本集，你會得到不靠譜的結果。我通常不會相信對任何少于50個樣本的問題的bootstrapping再抽樣，是以你最好也不要這樣做。
邊注，Seaborn的 barplot實際上使用bootstrapping來繪制置信區間：
seaborn.barplot(data=d, x='Month', y='Weight (kg)')
同樣，Seaborn非常适合進行探索性分析，其中一些圖表可以用于基本的統計。
回歸
讓我們提高一個檔次。我們将在這一群點上做一個線形回歸。
我将以我認為最普遍的方式來做這件事。我們将定義一個模型（在這種情況下是一條直線），一個損失函數（與該直線的平方偏差），然後使用通用求解器（scipy.optimize.minimize）對其進行優化。
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.plot(t_scale, model(x_scale, k_hat, m_hat), color='red',
linewidth=5, alpha=0.5)
具有不确定性的線性回歸，使用最大似然方法
我們隻拟合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)
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.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)
這裡的不确定性估計實際上不是100%，因為它沒有考慮k，m和σ本身的不确定性。這是一個不錯的近似，但要做到正确，我們需要同時做這些事情。是以我們這樣做吧。
再次啟用Bootstrapping！
是以讓我們把它提升到一個新的水準，并嘗試估計k和m和σ的不确定性估計！ 我認為這将顯示bootstrapping基本上是如何切割 - 你可以将它插入幾乎任何東西，以估計不确定性。
對于每個bootstrap估計，我将繪制一條線。我們也可以采用所有這些線并計算置信區間：
xys = list(zip(xs, ys))
curves = []
for i in range(100):
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)
哇，這裡發生了什麼？這種不确定性與之前的情節非常不同。這看起來很混亂，直到你意識到他們展示了兩個非常不同的東西：
● 第一個圖找到k和m的一個解，并顯示預測的不确定性。是以，如果你被問到下個月大象體重的範圍是什麼，你可以從圖表中得到它。
● 第二個圖找到了k和m的許多解，并顯示了kx + m的不确定性。是以，這回答了一個不同的問題 - 大象體重随時間變化的趨勢是什麼，趨勢的不确定性是什麼？
事實證明，我們可以将這兩種方法結合起來，并通過拟合繪制bootstrapping樣本并同時拟合k，m和σ使其更加複雜。然後，對于每個估算，我們可以預測新值y。我們這樣做：
for i in range(4000):
k_hat, m_hat, log_sigma_hat = scipy.optimize.minimize(
neg_log_likelihood, (0, 0, 0), args=(xs_bootstrap, ys_bootstrap)
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)
)
哎呦，又不錯！它現在變得很像回事了 - 如果仔細觀察，你會看到一個雙曲線形狀！
這裡的技巧是，對于（k，m，σ）的每個bootstrap估計，我們還需要繪制随機預測。正如您在代碼中看到的，我們實際上是将随機正态變量添加到y的預測值。這也是為什麼這個形狀最終變成一個大波浪形的原因。
不幸的是，bootstrapping對于這個問題來說相當緩慢 - 對于每個bootstrap，我們需要拟合一個模型。讓我們看看另一個選項：
馬爾可夫鍊蒙特卡羅方法
現在它會變得有點狂野了~~
我将切換到一些貝葉斯方法，我們通過繪制樣本來估計k，m和σ。它類似于bootstrapping，但MCMC有更好的理論基礎（我們使用貝葉斯規則從“後驗分布”中抽樣），并且它通常要快幾個數量級。
為此，我們将使用一個名為emcee的庫，我發現它很好用。 它所需要的隻是一個Log—likelihood函數，我們之前已經定義過了！ 我們隻需要用它的負數。
import emcee
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))
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')
你可以看到圍繞k = 200，m = 1000，σ= 100時的分布中心，我們原本也是如此建構它們的！這看起來挺令人放心的~
最後，我們可以使用與boostraps相同的方法繪制預測的完全不确定性：
samples = sampler.chain[:, -4000:, :].reshape((-1, ndim))
model(x_scale, k, m) +
numpy.exp(log_sigma) * numpy.random.normal(size=x_scale.shape)
這些貝葉斯方法并沒有在這裡結束。特别是有幾個庫可以解決這些問題。事實證明，如果您以更結構化的方式表達問題（而不僅僅是負對數似然函數），您可以将采樣比例調整為大問題（如數千個未知參數）。對于Python來說，有PyMC3和PyStan，以及稍微更小衆一點的Edward和Pyro。
結語
似乎我把你們帶進了一個深坑 - 但這些方法其實可以走得更遠。事實上，強迫自己估計我做的任何事情的不确定性是一個很好的體驗，可以逼着自己學習很多的統計知識。
根據資料做出決策很難！但如果我們對量化不确定性更加嚴格，我們可能會做出更好的決策。現在要做到這一點并不容易，但我真的希望我們能夠使用更便捷的工具，進而看到這些方法的普及。
原文釋出時間為：2018-11-27
本文來自雲栖社群合作夥伴“
大資料文摘”，了解相關資訊可以關注“
”。