貝葉斯A/B測試 - 一種計算兩種機率分佈差別性的方法過程

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()