import numpy as np
import pandas as pd

df1 = pd.DataFrame()
df1['treated'] = [94, 197, 16, 38, 99, 141, 23]
df2 = pd.DataFrame()
df2['control'] = [52, 104, 146, 10, 51, 30, 40, 27, 46]
df = pd.concat([df1, df2], axis=1)
df

print('estimated standard error by std / sqrt(n)')
print(df.std()/np.sqrt(df.count()))
print('estimated standard error by pandas function sem')
print(df.sem())  # standard error of the mean
se_diff = (df['treated'].sem()**2 + df['control'].sem()**2)**0.5
# for difference
m_diff = df['treated'].mean() - df['control'].mean()
print('difference of means = ', m_diff )
print('estimated standard error of the difference of means = ', se_diff)
print('difference / standard error of the difference = ', m_diff / se_diff)
df.describe()

import sys


tr = df['treated']
tr = tr[~np.isnan(tr)]  # NaN を省く.
cn = df['control']
size_tr = len(tr)
size_cn = len(cn)
print('size of tr and cn = ', size_tr, size_cn)
b = 2000
ms = []
ms_cn = []
for i in range(b):
    btr = np.random.choice(tr, size=size_tr)
    # print(btr)
    # sys.exit()
    m = btr.mean()
    ms.append(m)
    bcn = np.random.choice(cn, size=size_cn)
    m_cn = bcn.mean()
    ms_cn.append(m_cn)
ms = np.array(ms)
ms_cn = np.array(ms_cn)

dfm = pd.DataFrame()
dfm['treated'] = ms
dfm['control'] = ms_cn
dfm.hist(bins=20)

mms = ms.mean()
mms_cn = ms_cn.mean()

print('bootstrap estimate of the population means = ',  mms, mms_cn)

sqs_tr = np.sum( (ms - mms)**2 )
se_boot_tr = (sqs_tr / (b - 1))**.5
print(se_boot_tr)

sqs_cn = np.sum( (ms_cn - mms_cn)**2 )
se_boot_cn = (sqs_cn / (b - 1))**.5
print(se_boot_cn)

se_boot_difference = (se_boot_tr**2 + se_boot_cn**2)**.5
print(se_boot_difference)

print(' 処置  平均値    推定標準誤差')
print(' あり  {:.2f}   {:.2f}'.format(mms, se_boot_tr))
print(' なし  {:.2f}   {:.2f}'.format(mms_cn, se_boot_cn))
print(' 差　  {:.2f}   {:.2f}'.format(mms - mms_cn, se_boot_difference))
print(' 差/推定標準誤差  {:.2f}'.format((mms - mms_cn)/se_boot_difference))