Bayesian AB Testing Implemented in Two Methods 【Python】

# はじめに

この記事ではベイジアンABテストを2つの方法で実装しています。
1. Beta分布を共役事前分布に採用する方法 ([Adopt Beta distribution as the conjugate prior distribution](#adopt-beta-distribution-as-the-conjugate-prior-distribution))
2. Bernoulli分布のパラメータをベイズ推定する方法 ([Bayesian estimation of parameters of Bernoulli distribution](#bayesian-estimation-of-parameters-of-bernoulli-distribution))

## Check versions of libraries
```python
import numpy
import scipy
import pymc3



library_dict = {0: 'numpy', 1: 'scipy', 2: 'pymc3'}
for num, library in enumerate([numpy, scipy, pymc3]):
    print('{}: v{}'.format(library_dict[num], library.__version__))

>>
numpy: v1.18.5
scipy: v1.4.1
pymc3: v3.11.2
```

## Load libraries
```python
import numpy as np
from scipy.stats import beta, norm, bernoulli
import pymc3 as pm
import json
import sys
import logging
logger = logging.getLogger()
handler = logging.StreamHandler(sys.stdout)
handler.setLevel(logging.INFO)
logger.addHandler(handler)
logger.setLevel(logging.INFO)
from matplotlib import pyplot as plt
%matplotlib inline
from decimal import Decimal, ROUND_HALF_UP
```

## Adopt Beta distribution as the conjugate prior distribution
```python
class BayesianABN(object):
    
    def __init__(self, n_comparisons: int, params: list = None, n_samples: int = 50000) -> None:
        """
        n_comparisons: If u wanna compare two creatives, n_comparisons is 2.
        params: If u wanna compare two creatives, params is [[num of conversions, num of impressions], [num of conversions, num of impressions]].
        n_samples: Number of samples to generate from conjugate prior distribution.
        """
        if params is None:
            self.params = list()
            for _ in range(n_comparisons):
                self.params.append([1, 1])
        elif params is not None:
            self.params = params
            
        if n_comparisons != len(self.params):
            raise Exception("The value  of n_comparisons does not match the number of elements in params.")
        if n_comparisons * 2 != sum([len(list_) for list_ in self.params]):
            raise Exception("The value of n_comparisons does not match the number of elements contained in each element of params.")
            
        self.n_comparisons = n_comparisons
        self.data = [[0, 0] for _ in range(n_comparisons)]
        self.n_samples = n_samples
        self.sampling()
        
        logger.info(json.dumps({'Number of comparisons': self.n_comparisons, 'Number of beta samples': self.n_samples}))
        
    def sampling(self) -> None:
        self.posterior_distribution_list = [beta(*p).rvs(self.n_samples) for p in self.params]
        
    def update(self, additional_data) -> None:
        if self.n_comparisons * 2 != sum([len(list_) for list_ in additional_data]):
            raise Exception("The value of n_comparisons does not match the number of elements contained in each element of additional_data.")
        
        for num, data in enumerate(additional_data):
            cov = data[0]
            imp = data[1]
            
            self.data[num][0] += cov
            self.data[num][1] += imp
            
            self.params[num][0] += cov
            self.params[num][1] += imp - cov
            
        self.sampling()
        
    def significant_diff(self, target1: int, target2: int) -> float:
        self.diff = np.round((self.posterior_distribution_list[target1] <= self.posterior_distribution_list[target2]).mean(), 1)
        
        if self.diff < 0.5:
            return 1.0 - self.diff
            
        return self.diff
    
    def mean_var(self) -> dict:
        mean_var_dict = {chr(num + 65): {'mean': np.round(self.posterior_distribution_list[num].mean(), 2), 'var': Decimal(self.posterior_distribution_list[num].var()).quantize(Decimal('.000000'), rounding = ROUND_HALF_UP)} for num in range(self.n_comparisons)}
        
        return mean_var_dict
    
    def current_sample(self) -> dict:
        current_sample_dict = dict()
        for num in range(self.n_comparisons):
            current_sample_dict[chr(num + 65)] = dict()
            current_sample_dict[chr(num + 65)]['cv'] = self.data[num][0]
            current_sample_dict[chr(num + 65)]['imp'] = self.data[num][1]
        return current_sample_dict
        
    def distribution_drawing(self, title: str = '', save: bool = False, labels: str = None) -> None:
        plt.figure(figsize = (10, 5))
        plt.title(title)
        cmap = plt.get_cmap('jet')
        color_list = list()
        
        for num, posterior in enumerate(self.posterior_distribution_list):
            color = cmap(0.15 * (num + 1))
            color_list.append(color)
            plt.hist(posterior, bins = 100, histtype = 'stepfilled', density = True, color = color, alpha = 0.5)
        handles = [plt.Rectangle((0, 0), 1, 1, color = c, ec = 'k', alpha = 0.5) for c in color_list]
        
        if labels is None:
            labels = [chr(num + 65) for num in range(self.n_comparisons)]
            
        plt.legend(handles, labels)
        if save:
            plt.saving('{}.png'.format(title))
        plt.show()
        
    def delta_posterior_distribution_drawing(self, hdi_prob: float, target1: int, target2: int, title: str = None, save: bool = False) -> None:
        self.hdi_prob = hdi_prob
        self.pos_dis_1 = self.posterior_distribution_list[target1]
        self.pos_dis_2 = self.posterior_distribution_list[target2]
        
        self.delta_posterior_distribution = self.pos_dis_1 - self.pos_dis_2
        
        self.hdi = norm.interval(alpha = self.hdi_prob, loc = self.delta_posterior_distribution.mean(), scale = self.delta_posterior_distribution.std())
        
        plt.figure(figsize = (10, 5))
        if title is None:
            title = 'Difference in the posterior distributions of the parameters of {0} and {1} / {0} - {1} / {2}% HDI'.format(chr(target1 + 65), chr(target2 + 65), int(self.hdi_prob * 100))
        plt.title(title)
        cmap = plt.get_cmap('jet')
        color_list = list()
        
        color = cmap(0.15 * (0 + 1))
        color_list.append(color)
        plt.hist(self.delta_posterior_distribution, bins = 100, histtype = 'stepfilled', density = True, color = color, alpha = 0.5, label = 'delta')
        self.mean = self.delta_posterior_distribution.mean()
        plt.vlines(self.mean, 0, 10, linestyle = '-', label = 'mean: {}'.format(np.round(self.mean, 3)))
        if int(str(hdi_prob)[3]) == 0:
            plt.vlines(self.hdi[0], 0, 10, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round((1.0 - self.hdi_prob) / 2.0, 2) * 100, np.round(self.hdi[0], 3)))
            plt.vlines(self.hdi[1], 0, 10, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 2) * 100, np.round(self.hdi[1], 3)))
        elif int(str(hdi_prob)[3]) != 0:
            plt.vlines(self.hdi[0], 0, 10, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round((1.0 - self.hdi_prob) / 2.0, 3) * 100, np.round(self.hdi[0], 3)))
            plt.vlines(self.hdi[1], 0, 10, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 3) * 100, np.round(self.hdi[1], 3)))
        handles = [plt.Rectangle((0, 0), 1, 1, color = c, ec = 'k', alpha = 0.5) for c in color_list]
            
        plt.legend()
        if save:
            plt.saving('{}.png'.format(title))
        plt.show()
```
```python
BayesianABN_instance = BayesianABN(3)

>> {"Number of comparisons": 3, "Number of beta samples": 50000}

BayesianABN_instance.update([[2, 204], [3, 200], [50, 300]])

BayesianABN_instance.significant_diff(0, 1)

>> 0.7

BayesianABN_instance.current_sample()

>>
{'A': {'cv': 2, 'imp': 204},
 'B': {'cv': 3, 'imp': 200},
 'C': {'cv': 50, 'imp': 300}}

BayesianABN_instance.mean_var()

>> 

{'A': {'mean': 0.01, 'var': Decimal('0.000070')},
 'B': {'mean': 0.02, 'var': Decimal('0.000094')},
 'C': {'mean': 0.17, 'var': Decimal('0.000461')}}
```
```pytyon
BayesianABN_instance.distribution_drawing()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/b2892258-20d2-b74e-701f-3bc3bbc579ff.png">

## Bayesian estimation of parameters of Bernoulli distribution
```python
class BayesianAB(object):
    
    def __init__(self, sample_dict: dict):
        """
        sample_dict: If u wanna compare two creatives, sample_dict = {'A': creative_a, 'B': creative_b}. The values in the dict contain a sampling of cv = 1/ not cv = 0 for each creative.
        """
        self.sample_dict = sample_dict
        if len(self.sample_dict) != 2:
            raise Exception("The number of elements in sample_dict must be two.")
    
    def modeling(self) -> None:
        with pm.Model() as self.model:
            self.p_A = pm.Uniform('p_A', 0, 1.0)
            self.p_B = pm.Uniform('p_B', 0, 1.0)
            
            self.obs_A = pm.Bernoulli('obs_A', self.p_A, observed = self.sample_dict['A'])
            self.obs_B = pm.Bernoulli('obs_B', self.p_B, observed = self.sample_dict['B'])
            
            self.delta_prob = pm.Deterministic('delta_prob', self.p_B - self.p_A)
            
    def estimate(self) -> None:
        with self.model:
            self.start = pm.find_MAP()
            self.step = pm.Slice()
            self.trace = pm.sample(20000, step = self.step, start = self.start, return_inferencedata = False)
            self.burned_trace = self.trace[1000:]
            
    def model_to_graphviz(self) -> None:
        return pm.model_to_graphviz(self.model)
    
    def tracing(self) -> None:
        with self.model:
            pm.plot_trace(self.trace)
    
    def posterior_distribution_drawing(self, hdi_prob) -> None:
        with self.model:
            pm.plot_posterior(self.burned_trace, hdi_prob = hdi_prob)
        
    def summary(self, hdi_prob: float):
        with self.model:
            self.summary_ = pm.summary(self.burned_trace, hdi_prob = hdi_prob)
            return self.summary_
        
    def significant_diff(self) -> float:
        self.diff = np.round((self.burned_trace['p_A'] <= self.burned_trace['p_B']).mean(), 1)
        
        if self.diff < 0.5:
            return 1.0 - self.diff
            
        return self.diff
    
    def delta_posterior_distribution_drawing(self, hdi_prob: float, title: str = None, save: bool = False) -> None:
        self.hdi_prob = hdi_prob
        try:
            self.summary_
        except:
            with self.model:
                self.summary_ = pm.summary(self.burned_trace, hdi_prob = self.hdi_prob)
        delta_prob = self.burned_trace['delta_prob']
        
        plt.figure(figsize = (10, 5))
        if title is None:
            title = 'Posterior distribution of delta'
        plt.title(title)
        
        cmap = plt.get_cmap('jet')
        color_list = list()
        color = cmap(0.15 * (0 + 1))
        color_list.append(color)
        
        plt.vlines(self.summary_.iloc[2]['mean'], 0, 50, linestyle = '-', label = 'mean: {}'.format(self.summary_.iloc[2]['mean']))
        if int(str(hdi_prob)[3]) == 0:
            self.hdi_lower = self.summary_.iloc[2]['hdi_{}%'.format(np.round((1.0 - self.hdi_prob) / 2.0, 2) * 100)]
            self.hdi_upper = self.summary_.iloc[2]['hdi_{}%'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 2) * 100)]
            plt.vlines(self.hdi_lower, 0, 50, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round((1.0 - self.hdi_prob) / 2.0, 2) * 100, self.hdi_lower))
            plt.vlines(self.hdi_upper, 0, 50, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 2) * 100, self.hdi_upper))
        elif int(str(hdi_prob)[3]) != 0:
            self.hdi_lower = self.summary_.iloc[2]['hdi_{}%'.format(np.round((1.0 - self.hdi_prob) / 2.0, 3) * 100)]
            self.hdi_upper = self.summary_.iloc[2]['hdi_{}%'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 3) * 100)]
            plt.vlines(self.hdi_lower, 0, 50, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round((1.0 - self.hdi_prob) / 2.0, 3) * 100, self.hdi_lower))
            plt.vlines(self.hdi_upper, 0, 50, linestyle = '--', label = 'hdi_{}%: {}'.format(np.round(1.0 - (1.0 - self.hdi_prob) / 2.0, 3) * 100, self.hdi_upper))
        plt.hist(delta_prob, bins = 100, histtype = 'stepfilled', density = True, color = color, alpha = 0.5)
        handles = [plt.Rectangle((0, 0), 1, 1, color = c, ec = 'k', alpha = 0.5) for c in color_list]
            
        plt.legend()
        if save:
            plt.saving('{}.png'.format(title))
        plt.show()
```

## [Adopt Beta distribution as the conjugate prior distribution](#adopt-beta-distribution-as-the-conjugate-prior-distribution) VS [Bayesian estimation of parameters of Bernoulli distribution](#bayesian-estimation-of-parameters-of-bernoulli-distribution)
### Create Dummy Data
```python
creative_a = bernoulli.rvs(p = 0.15, size = 15000)
creative_b = bernoulli.rvs(p = 0.17, size = 13000)
```

### Comparison
#### Adopt Beta distribution as the conjugate prior distribution
```python
BayesianABN_instance = BayesianABN(2)

>> {"Number of comparisons": 2, "Number of beta samples": 50000}

BayesianABN_instance.update([[sum(creative_a), len(creative_a)], [sum(creative_b), len(creative_b)]])
BayesianABN_instance.significant_diff(0, 1)

>> 1.0

BayesianABN_instance.current_sample()

>> {'A': {'cv': 2258, 'imp': 15000}, 'B': {'cv': 2244, 'imp': 13000}}
```
```python
BayesianABN_instance.distribution_drawing()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/545fb098-ec84-4436-d246-bddfe3383743.png">

```python
BayesianABN_instance.delta_posterior_distribution_drawing(0.95, 0, 1)
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/b8dd6355-6ee3-4581-3363-fe9bd738959b.png">

#### Bayesian estimation of parameters of Bernoulli distribution
```python
instance = BayesianAB({'A': creative_a, 'B': creative_b})
instance.modeling()
instance.model

>> 
p_A_interval__ ~ TransformedDistribution
p_B_interval__ ~ TransformedDistribution
           p_A ~ Uniform
           p_B ~ Uniform
    delta_prob ~ Deterministic
         obs_A ~ Bernoulli
         obs_B ~ Bernoulli
```
```python
instance.model_to_graphviz()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/b8db044f-ccaf-dbc3-274d-a46518d54404.png">

```python
instance.estimate()

>> 
 100.00% [7/7 00:00<00:00 logp = -19,408, ||grad|| = 6,752.2]

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>Slice: [p_B]
>Slice: [p_A]

 100.00% [84000/84000 01:15<00:00 Sampling 4 chains, 0 divergences]
Sampling 4 chains for 1_000 tune and 20_000 draw iterations (4_000 + 80_000 draws total) took 120 seconds.
```
```python
instance.tracing()
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/890c007c-b459-d03c-64ca-042520739261.png">

```python
posterior_distributions = instance.posterior_distribution_drawing(hdi_prob = 0.95)
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/7532a628-9923-00df-5d06-b0cbccb3fe39.png">

```python
summary = instance.summary(hdi_prob = 0.95)
summary

>> 
           mean   sd   hdi_2.5%   hdi_97.5%   mcse_mean   mcse_sd   ess_bulk   ess_tail   r_hat
p_A        0.151 0.003   0.145      0.156        0.0        0.0      76771.0   58106.0     1.0
p_B        0.173 0.003   0.166      0.179        0.0        0.0      73469.0   56450.0     1.0
delta_prob 0.022 0.004   0.013      0.031        0.0        0.0      75624.0   66133.0     1.0
```
```python
instance.significant_diff()

>> 1.0

(instance.burned_trace['p_B'] - instance.burned_trace['p_A'] > 0).mean()

>> 1.0
```
```python
instance.delta_posterior_distribution_drawing(hdi_prob = 0.95, title = 'Posterior distribution of delta / 95% HDI')
```
<img src = "https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/696827/6f399a40-1c1e-2e28-c778-150f88ba2d09.png">

#### Result
The same average, HID2.5%, and HDI 97.5% were calculated for the two methods compared. 
Since the difference of the estimated parameters is 1.0 and the difference of the average is 0.022, we can say that the probability that B is more likely to convert is close to 100% and the conversion rate could have been increased by 2.2%.

# まとめ
ベイジアンABテストを2つの方法で実装し結果を比較してみました。
認識や実装の間違いなどありましたらご指摘のほどよろしくお願いいたします。