1. 控制變數
0x1：控制變數主要思想
科學中對於多因素（多變數）的問題，常常採用控制因素（變數）的方法，吧多因素的問題變成多個單因素的問題。每一次只改變其中的某一個因素，而控制其餘幾個因素不變，從而研究被改變的這個因素對事物的影響，分別加以研究，最後再綜合解決，這種方法叫控制變數法。
它是科學探索中的重要思想方法，廣泛地運用在各個科學探索和科學實驗研究之中。
0x2：控制變數思想在機器學習中的應用
在機器學習專案中，我們可能會將專家領域經驗融合到特徵工程中，即主觀先驗。
在設計並獲得特徵向量後，我們會在直方圖上列印出pdf概率密度圖，目的是檢視不同特徵之間的區分度。
這個比較過程，本質上就是控制變數過程，在該實驗中，自變數是兩個不同的特徵，因變數是目標值（label值），其他因素都完全相同，通過對比兩個不同特徵的概率分佈差異性（例如使用t-test），來得到特徵可區分型的判斷。
如果兩個特徵之間差異性很低，則說明可能需要進行feature selection後者feature reduction。
0x3：控制變數思想在貝葉斯統計推斷中的作用
在A/B測試中，兩組實驗分別採用了不同的策略，其他因素完全相同，通過對實驗的後驗結果進行差異性分佈（例如t-test），以此得到這兩種策略之間是否存在明顯的差異或者說增益的判斷（判斷本身也是概率性的）。
Relevant Link:
https://wenku.baidu.com/view/afe92b8fb04e852458fb770bf78a6529647d35c7.html
2. A/B測試簡介
0x1：A/B測試基本思想
A/B測試背後的基本思想是：假如有一個理想的平行宇宙用於對照，該宇宙的人與我們這裡的人是完全一樣的（除控制變數外，其他變數都相同），那麼此時如果給某一邊的人以某種特殊的待遇（調整控制變數），那麼結果所導致的變化一定會被歸咎於這一特殊待遇。
但是在實踐中，我們沒法進入到平行宇宙，因此我們只能利用兩組足夠大量的樣本來近似地創造一對平行宇宙。
0x2：A/B測試的最終目的
A/B測試的最終目的是比較A/B的後驗分佈的差異性，使用的方法有很多種，例如t-test。
0x3：A/B測試的優勢
1. 可解釋的概率值。在貝葉斯分佈中，我們可以直接回答諸如”我們出錯的概率是多少“之類的問題，這在頻率派的方法中通常難以回答； 2. 很容易應用損失函式；
我們在這篇blog中用一個網站轉化率測試的簡單案例，來一起討論下A/B測試。
Relevant Link:
https://baike.baidu.com/item/AB%E6%B5%8B%E8%AF%95/9231223?fr=aladdin
3. 二項分佈問題場景的A/B測試 - 一個網站轉化率測試案例
我們有A和B兩種網站設計。當使用者登入網站時，我們隨機地將其引向其中之一（模擬隨機分組），並且記錄下來。當有足夠多的使用者訪問以後，我們用得到的資料來計算兩種設計的轉化率指標。
考慮我們得到了如下數字：
visitors_to_A = 1300； visitors_to_B = 1300； conversions_from_A = 120； conversions_from_B = 125；
我們關心的是A和B的轉化概率。從商業化角度考慮，我們希望轉化率越高越好。
為此，我們需要對A和B的轉化率進行建模。
0x1：選擇先驗分佈
由於需要對概率建模，因此選擇Beta分佈作為先驗是一個好主意，轉化率的取值範圍在0~1之間。
同時，我們的訪客數量和轉化資料是二項式的，每個訪客只有兩種可能：轉化 or 不轉化。
Beta分佈的共軛特性是的我們不需要進行MCMC，可以直接得到後驗概率分佈。
0x2：計算後驗分佈
假設我們的先驗是Beta（1，1），它等價於【0，1】上的均勻分佈。
輸入樣本資料後，得到Beta後驗分佈：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from scipy.stats import beta visitors_to_A = 1300 visitors_to_B = 1300 conversions_from_A = 120 conversions_from_B = 125 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversions_from_A, beta_prior + visitors_to_A - conversions_from_A) posterior_B = beta(alpha_prior + conversions_from_B, beta_prior + visitors_to_B - conversions_from_B) plt.show()
0x3：比較後驗分佈的大小
1. MCMC後驗取樣法
接下來我們想判斷哪個組轉化率可能更高。為此，類似MCMC的做法，我們對後驗進行取樣，並且比較來自於A的後驗樣本的概率大於來自B的後驗樣本的概率：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from scipy.stats import beta visitors_to_A = 1300 visitors_to_B = 1300 conversions_from_A = 120 conversions_from_B = 125 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversions_from_A, beta_prior + visitors_to_A - conversions_from_A) posterior_B = beta(alpha_prior + conversions_from_B, beta_prior + visitors_to_B - conversions_from_B) samples = 2000 samples_posterior_A = posterior_A.rvs(samples) samples_posterior_B = posterior_B.rvs(samples) print(samples_posterior_A > samples_posterior_B).mean()
可以看到，有31%的概率A比B的轉化效率高。或者說，有69%的概率B比A的轉化率高。
這並不是十分顯著，因為如果兩個頁面完全相同，那麼重新實驗得到的概率會接近50%，而我們實驗中得出的結論比較接近50%，這個結論並不能提供多少有用的資訊。
通過概率密度圖顯示兩個分佈的結果：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from scipy.stats import beta from IPython.core.pylabtools import figsize visitors_to_A = 1300 visitors_to_B = 1300 conversions_from_A = 120 conversions_from_B = 125 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversions_from_A, beta_prior + visitors_to_A - conversions_from_A) posterior_B = beta(alpha_prior + conversions_from_B, beta_prior + visitors_to_B - conversions_from_B) samples = 2000 samples_posterior_A = posterior_A.rvs(samples) samples_posterior_B = posterior_B.rvs(samples) print(samples_posterior_A > samples_posterior_B).mean() figsize(12.5, 4) plt.rcParams['savefig.dpi'] = 300 plt.rcParams['figure.dpi'] = 300 x = np.linspace(0, 1, 500) plt.plot(x, posterior_A.pdf(x), label='posterior of A') plt.plot(x, posterior_B.pdf(x), label='posterior of B') plt.xlim(0.05, 0.15) plt.xlabel('value') plt.ylabel('density') plt.show()
2. 通過KL散度計算後驗分佈的差異性
Relevant Link:
https://baike.baidu.com/item/%E7%9B%B8%E5%AF%B9%E7%86%B5/4233536?fromtitle=KL%E6%95%A3%E5%BA%A6&fromid=23238109&fr=aladdin
4. 多項分佈下的A/B測試
網際網路公司的一個常見的目標是，不僅是要增加註冊量，還要優化使用者可能選擇的註冊方案。比如，一個業務體可能希望使用者在多種備選項中，選擇價格更高的方案。
假設給使用者展示兩個版本的定價頁面，並且我們希望得到每次訪問的收入期望值，不僅僅是關心使用者是否註冊，而是想知道能獲得的收入的期望值。
0x1：收入期望的分析
企業網頁總收入的期望為：
這裡，P79 為選擇 $79 收費方案的概率，其他的類似，$0 代表使用者未選擇任何收費方案（既為轉化）。這樣一來，整體概率和為1：
0x2：選擇概率分佈模型
接下來是為各個收費方案的選擇先驗分佈。這裡不能簡單使用Beta分佈，因為各個選項之間彼此是相關的（不符合beta分佈成立條件），它們的和為1。比如，如果 P79 很大，那麼其他的概率必然較小。
Beta分佈有一個推廣，叫做 Dirichlet（狄利克雷）分佈。它返回一個和為 1 的正數陣列。陣列的長度由一個輸入向量的長度決定，這一輸入向量的值類似於先驗的引數。
另一方面，對於樣本資料的建模，也不能選擇二項分佈，因為選項不只2個。二項分佈有一個推廣叫做多項分佈，我們的觀測值服從多項分佈，並且各個取值的概率都是未知的。
同時，狄利克雷分佈是多項分佈的共軛先驗！這意味著對於位置概率的後驗，我們有明確的公式。
如果先驗服從 Dirichlet（1，1，...，1），並且我們的觀測值為 N1，N2，....，Nm，那麼後驗是：
0x3：計算後驗分佈
假如有1000人瀏覽了頁面，並且註冊情況如下：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from numpy.random import dirichlet from IPython.core.pylabtools import figsize N = 1000 N_79 = 10 N_49 = 46 N_25 = 80 N_0 = N - (N_79 + N_49 + N_25) observations = np.array([N_79, N_49, N_25, N_0]) prior_parameters = np.array([1, 1, 1, 1]) posterior_samples = dirichlet(prior_parameters+observations, size=10000) print "Two random samples from the posterior: " print posterior_samples[0] print posterior_samples[1] for i, label in enumerate(['p_79', 'p_49', 'p_25', 'p_0']): ax = plt.hist(posterior_samples[:, i], bins=50, label=label, histtype='stepfilled') plt.xlabel('value') plt.ylabel('density') plt.show()
從上圖可以看出，關於概率的可能取值仍然有不確定性，所有期望值的結果也是不確定的。我們得到的就是期望值的後驗分佈。
0x4：計算總收入期望的後驗分佈
來自該後驗的樣本的和總是1，因此可以將這些樣本用於前面的期望值的公式裡參與計算收入的後驗期望。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from numpy.random import dirichlet from IPython.core.pylabtools import figsize N = 1000 N_79 = 10 N_49 = 46 N_25 = 80 N_0 = N - (N_79 + N_49 + N_25) observations = np.array([N_79, N_49, N_25, N_0]) prior_parameters = np.array([1, 1, 1, 1]) posterior_samples = dirichlet(prior_parameters+observations, size=10000) def expected_revenue(P): return 79 * P[:, 0] + 49 * P[:, 1] + 25 * P[:, 2] + 0 * P[:, 3] posterior_expected_revenue = expected_revenue(posterior_samples) plt.hist(posterior_expected_revenue, histtype='stepfilled', label='expected revenue', bins=50) plt.xlabel('value') plt.ylabel('density') plt.show()
從上圖中我們可以解讀出幾個資訊：
1. 收入的期望值有很大可能在 $4 和 $6 之間，不大可能在這個範圍之外
0x5：對兩個不同的頁面進行分析 - 延伸到A/B測試
接下來我們把問題擴充套件一些，嘗試對兩個不同的WEB頁面進行這樣的分佈，從而進行A/B測試，對比在相同的條件下，A、B頁面的後驗分佈的概率分佈差異性。
1. 計算後驗分佈
我們將兩個站點稱為A和B，併為它們虛構一些資料：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from numpy.random import dirichlet from IPython.core.pylabtools import figsize N_A = 1000 N_A_79 = 10 N_A_49 = 46 N_A_25 = 80 N_A_0 = N_A - (N_A_79 + N_A_49 + N_A_25) observations_A = np.array([N_A_79, N_A_49, N_A_25, N_A_0]) N_B = 2000 N_B_79 = 45 N_B_49 = 84 N_B_25 = 200 N_B_0 = N_B - (N_B_79 + N_B_49 + N_B_25) observations_B = np.array([N_B_79, N_B_49, N_B_25, N_B_0]) prior_parameters = np.array([1, 1, 1, 1]) posterior_samples_A = dirichlet(prior_parameters+observations_A, size=10000) posterior_samples_B = dirichlet(prior_parameters+observations_B, size=10000) def expected_revenue(P): return 79 * P[:, 0] + 49 * P[:, 1] + 25 * P[:, 2] + 0 * P[:, 3] posterior_expected_revenue_A = expected_revenue(posterior_samples_A) posterior_expected_revenue_B = expected_revenue(posterior_samples_B) plt.hist(posterior_expected_revenue_A, histtype='stepfilled', label='expected revenue A', bins=50) plt.hist(posterior_expected_revenue_B, histtype='stepfilled', label='expected revenue B', bins=50, alpha=0.8) plt.xlabel('value') plt.ylabel('density') plt.show()
2. 從後驗分佈概率圖中我們得到了什麼初步的推斷資訊
從上圖中我們可以得到如下資訊：
1. 兩個後驗的距離比較遠，說明兩種頁面的表現有很多差別，但這裡只是一個感性的認知，具體差別多少不知道； 2. 頁面A的累計期望比頁面B要少 1$，看起來不多，但是每次瀏覽都有這樣的差距，積少成多是很可觀的，越是大型企業對這點感受越深；
3. 定量分析兩種頁面的後驗差距
為了確認這種差距是真實存在的，我們來看看看看概率密度差距：
print(posterior_expected_revenue_B > posterior_expected_revenue_A).mean() 0.9804
結果為98%，這個值已經足夠高了，業務方應該選擇頁面B的方案。
另一個有趣的視角是兩種頁面的後驗差距，我們需要看看兩種收入期望的後驗在直方圖中的間距：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from numpy.random import dirichlet from IPython.core.pylabtools import figsize N_A = 1000 N_A_79 = 10 N_A_49 = 46 N_A_25 = 80 N_A_0 = N_A - (N_A_79 + N_A_49 + N_A_25) observations_A = np.array([N_A_79, N_A_49, N_A_25, N_A_0]) N_B = 2000 N_B_79 = 45 N_B_49 = 84 N_B_25 = 200 N_B_0 = N_B - (N_B_79 + N_B_49 + N_B_25) observations_B = np.array([N_B_79, N_B_49, N_B_25, N_B_0]) prior_parameters = np.array([1, 1, 1, 1]) posterior_samples_A = dirichlet(prior_parameters+observations_A, size=10000) posterior_samples_B = dirichlet(prior_parameters+observations_B, size=10000) def expected_revenue(P): return 79 * P[:, 0] + 49 * P[:, 1] + 25 * P[:, 2] + 0 * P[:, 3] posterior_expected_revenue_A = expected_revenue(posterior_samples_A) posterior_expected_revenue_B = expected_revenue(posterior_samples_B) posterior_diff = posterior_expected_revenue_B - posterior_expected_revenue_A plt.hist(posterior_diff, histtype='stepfilled', color='#7A68A6', label='difference in revenue between B and A', bins=50) plt.vlines(0, 0, 700, linestyles='--') plt.xlabel('value') plt.ylabel('density') plt.show()
從上圖中可以看到：
1. 兩者兼具有50%的概率大於 1$，並且有一定可能大於 2$； 2. 即使我們選擇B是錯誤的（虛線左邊），也不會有太大的損失，分佈上幾乎不會超出 -0.5$ 太多；
0x6：A/B兩組後驗的增幅的估計
在進行了A/B測試後，決策者通常會對增幅感興趣。但實際上，這裡用增幅這個詞是不準確的，貝葉斯估計的後驗結果是一個分佈，兩個分佈的增幅也是一個分佈。
我們在學習貝葉斯統計的時候，一定要時刻注意將連續值問題和二值問題區分開來：
1. 連續值問題是要衡量結果到底好多少，這是一定範圍內的連續值（軟分界）； 2. 二值問題是要判斷誰更好，只有兩種可能（硬分界）；
1. 用後驗分佈的均值計算相對增幅合理嗎？
一個很自然的想法是，用兩個後驗的均值計算相對增幅：
這會帶來一些嚴重的錯誤。首先，這把對 Pa 和 Pb 的真實值的不確定性都掩蓋起來了。在用均值差公式來計算增幅時，我們假定了這些值都是精確已知的（均值本質是一個統計壓縮方法）。這幾乎總是會高度這些值。
2. 計算後驗分佈的概率密度增幅
解決上述問題的方法就是，保留不確定性，統計學畢竟就是關於不確定性的理論。為此，我們可以直接將後驗分佈的概率密度函式相減，得到一個相減後的概率分佈。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from scipy.stats import beta from IPython.core.pylabtools import figsize visitors_to_A = 1300 visitors_to_B = 1300 conversions_from_A = 120 conversions_from_B = 125 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversions_from_A, beta_prior + visitors_to_A - conversions_from_A) posterior_B = beta(alpha_prior + conversions_from_B, beta_prior + visitors_to_B - conversions_from_B) samples = 2000 samples_posterior_A = posterior_A.rvs(samples) samples_posterior_B = posterior_B.rvs(samples) print(samples_posterior_A > samples_posterior_B).mean() def relative_increase(a, b): return (a - b) / b posterior_rel_increase = relative_increase(samples_posterior_A, samples_posterior_B) figsize(12.5, 4) plt.hist(posterior_rel_increase, label='relative increase') plt.xlabel('value') plt.ylabel('density') plt.show()
從上圖中可以看出：
1. 有89%的可能性，相對增幅會達到20%或者更多； 2. 有72%的可能性，增幅能達到50%；
如果我們想要簡單地使用均值點估計，即：
則關於增幅的均值點估計應該是87%，這顯然太高了。
3. 建立後驗概率增幅的點估計
我們繼續增幅計算這個話題的討論，儘管從貝葉斯統計角度來說，一切都是概率。但是直接把一個分佈交給業務方是不合適的，業務方希望得到的結果就是一個精確的數值。那怎麼辦呢？解決的方法就是用統計壓縮的方式，從分佈中提取出一個”有代表性“的精確統計量來代替分佈，有三種可選的方案：
1. 返回增幅後驗分佈的均值：這個方法並不是非常好。原因在於對於一個傾斜的長尾分佈，類似均值這樣的統計量會很受影響，因而結論會過分表達長尾資料以至於高估實際的相對增幅； 2. 返回增幅後驗分佈的中位數：中位數應該是更合理的值，它對於傾斜的分佈會更有魯棒性。然而在實踐中，中位數仍然可能導致結論被高估； 3. 返回增幅後驗分佈的百分位數（低於50%）。比如，返回第30%百分位數。這樣做會有兩個想要的特性： 1）它相當於從數學上在增幅後驗分佈之上應用了一個損失函式。以懲罰過高的估計，這樣估計的結果就更加保守 2）隨著我們得到越來越多的實驗資料，增幅的後驗分佈會越來越窄，意味著任何百分位數都會收斂到同一個點
我們在下圖中把三種統計量都畫出來：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import pymc as pm from scipy.stats import beta from IPython.core.pylabtools import figsize visitors_to_A = 1275 visitors_to_B = 1300 conversions_from_A = 22 conversions_from_B = 12 alpha_prior = 1 beta_prior = 1 posterior_A = beta(alpha_prior + conversions_from_A, beta_prior + visitors_to_A - conversions_from_A) posterior_B = beta(alpha_prior + conversions_from_B, beta_prior + visitors_to_B - conversions_from_B) samples = 20000 samples_posterior_A = posterior_A.rvs(samples) samples_posterior_B = posterior_B.rvs(samples) print(samples_posterior_A > samples_posterior_B).mean() def relative_increase(a, b): return (a - b) / b posterior_rel_increase = relative_increase(samples_posterior_A, samples_posterior_B) mean = posterior_rel_increase.mean() median = np.percentile(posterior_rel_increase, 50) conservtive_percentile = np.percentile(posterior_rel_increase, 30) figsize(12, 4) plt.hist(posterior_rel_increase, label='relative increase') plt.vlines(mean, 0, 2500, linestyles='-.', label='mean') plt.vlines(median, 0, 2500, linestyles=':', label='median') plt.vlines(conservtive_percentile, 0, 2500, linestyles='--', label='conservtive_percentile') plt.xlabel('value') plt.ylabel('density') plt.show()
