我花了一年時間研究不確定性估算，寫下了這份最全指南
讓我從年初的一個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 group s:
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()
神奇吧，這些圖表與之前的圖表非常相似！（正如我們本該期待的那樣）
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)' )
哇，這裡發生了什麼？這種不確定性與之前的情節非常不同。這看起來很混亂，直到你意識到他們展示了兩個非常不同的東西：
● 第一個圖找到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。
結語
似乎我把你們帶進了一個深坑 - 但這些方法其實可以走得更遠。事實上，強迫自己估計我做的任何事情的不確定性是一個很好的體驗，可以逼著自己學習很多的統計知識。
根據資料做出決策很難！但如果我們對量化不確定性更加嚴格，我們可能會做出更好的決策。現在要做到這一點並不容易，但我真的希望我們能夠使用更便捷的工具，從而看到這些方法的普及。
原文釋出時間為：2018-11-27