二項分布に従う標本から施策の効果(ATE)を推定

仕事で始める機械学習７章より。一部改変。

```python
import numpy as np
import pandas as pd
import japanize_matplotlib
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import scipy.stats
import scipy
import statsmodels.tools.print_version
```

## ATE=0.015として二項分布に従う標本を抽出。
施策の効果を推定量の差と対照群に対する効果(対照群での成約率のこと)を推定量の従う分布とその信頼区間から算出し、施策継続の可否を判断。

```python
# 二項分布に従う対照群と介入群の標本を抽出
T = 10000
# Random variates of given type.
random_state=0

control = scipy.stats.binom(1, 0.05).rvs(size=T,random_state=random_state)
treatment = scipy.stats.binom(1, 0.065).rvs(size=T,random_state=random_state)
print('E(Y0)=', np.mean(control),',','E(Y1)=', np.mean(treatment))

n_observed = np.ones(T).cumsum().astype(int)
treatment_p = treatment.cumsum()/n_observed
control_p = control.cumsum()/n_observed

```

    E(Y0)= 0.0494 , E(Y1)= 0.0642
    


```python
import statsmodels.stats.proportion
# @np.vectorize
def calc_confidence_interval(n_success, n_observed):
   return statsmodels.stats.proportion.proportion_confint(
      n_success, n_observed, alpha=0.05, method='wilson'
   )

treatment_lcb, treatment_ucb = calc_confidence_interval(treatment.cumsum(), n_observed)
control_lcb, control_ucb = calc_confidence_interval(control.cumsum(), n_observed)

fig = plt.figure(figsize=(8, 5))
plt.plot(n_observed, treatment_p, label='介入群の推定量')
plt.fill_between(n_observed, treatment_lcb, treatment_ucb, alpha=0.3, interpolate=True)
plt.plot(n_observed, control_p, linestyle='dashed',  label='対照群の推定量')
plt.fill_between(n_observed, control_lcb, control_ucb, alpha=0.3, interpolate=True)
plt.legend(loc='best')
plt.xlim(0, 10000)
plt.xlabel('観測データの数')
plt.ylabel('E(Y)')
plt.ylim(0.030, 0.08)
```




    (0.03, 0.08)




    
![output_3_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1347377/5675c224-3191-53ad-a6e1-cb0d071b12b6.png)
    


- 2群の差の信頼区間の下限が0より高ければATEあり、求める施策効果の水準によって施策の採択を決定する
- 例えば対照群と比較して0.13~0.38%の有料契約率の向上が見込める例で、月契約ページ閲覧数100000人、月額料金980円(税抜き)とすると
    - 100000人に対して5000人契約が5650~6900人であるので、
    - 650*980=637000 ~ 1900*980=1862000
- よって施策の効果が637000~1862000であるとわかる。




```python
from statsmodels.stats.weightstats import DescrStatsW, CompareMeans

def _calc_diff_confidence_interval(n_success_a, n_success_b, n_observed):
    if n_observed < 10:
        return (1, -1)
    # print(n_success_a)
    # print(np.ones(n_success_a))
    # print(n_observed - n_success_a)
    # print(np.zeros(n_observed - n_success_a))
    # print(np.append(np.ones(n_success_a), np.zeros(n_observed - n_success_a)))
    sample_a = DescrStatsW(np.append(np.ones(n_success_a), np.zeros(n_observed - n_success_a)))
    sample_b = DescrStatsW(np.append(np.ones(n_success_b), np.zeros(n_observed - n_success_b)))
    # print(CompareMeans(sample_a, sample_b).zconfint_diff(alpha=0.05))
    lcb, ucb = CompareMeans(sample_a, sample_b).zconfint_diff(alpha=0.05)
    
    return float(lcb), float(ucb)

calc_diff_confidence_interval = np.vectorize(_calc_diff_confidence_interval, otypes=[np.float, np.float])

diff_interval_lcb, diff_interval_ucb = calc_diff_confidence_interval(
   treatment.cumsum(),
   control.cumsum(),
   n_observed
)

diff = treatment_p - control_p

xmax = 10000
fig = plt.figure(figsize=(8, 5))
plt.plot(n_observed, diff, color='green', label='介入効果の推定量')
plt.fill_between(n_observed, diff_interval_lcb, diff_interval_ucb, alpha=0.6, color='lightgreen', interpolate=True)
plt.hlines([0], 0, xmax, "black")
plt.hlines([0.015], 0, xmax, "green", linestyles='dashed', label='真の介入効果')
plt.ylim(-0.02, 0.03)
plt.ylabel('2群のE(Y)の差')
plt.xlabel('観測データの数')
plt.legend(loc='best')
plt.xlim(0, xmax)
```




    (0.0, 10000.0)




    
![output_5_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1347377/6c515bc2-1e6c-b981-1565-5710277cb319.png)
    



```python
print('対照群との比較で申し込み率の向上は{:.2g}~{:.2g}%です。'.format(diff_interval_lcb[-1]/np.mean(control),diff_interval_ucb[-1]/np.mean(control)))
```
    対照群との比較で申し込み率の向上は0.17~0.43%です。