Return to Revisions

Bootstrap

I tried a bootstrap. It often threw the error:

division by zero

I got around that by assuming:

infinity = approximately 10101

with the resulting estimates:

Figure 1, Improvement (I) against Bootstrap Size (BS) for runs = 100

Figure 2, Improvement (I) against against Runs for BS=28

This gave an improvement after 1000 runs with BS=28 of:

55% [14%, 150%]

In Python:

import numpy  as np
import pandas as pd
import pickle
import matplotlib.pyplot as plt

print('Generate the sample data')
data = pd.DataFrame({'A':[1]*11+[0]*6+[0]*11,
                     'B':[1]*11+[1]*6+[0]*11})
print('sample size: ',len(data))
print('')
print('A B X')
print('1 1',len(data[((data.A==1)&(data.B==1))]))
print('1 0',len(data[((data.A==1)&(data.B==0))]))
print('0 1',len(data[((data.A==0)&(data.B==1))]))
print('0 0',len(data[((data.A==0)&(data.B==0))]))
print('')

# Results
I     = {}
Lower = {}
Media = {}
Upper = {}

# Control Parameters
Runs = range(100)
#bootstrap_size = range(len(data))
BS_Max = 100
bootstrap_size = range(BS_Max)

for BS in bootstrap_size:
    #print('bootstrap size ', BS)
    
    # Results
    I_T = {}
    
    for R in Runs:
        
        # Bootstrap
        BooP = data.sample(BS, replace=True)
        
        # Data
        X_11 = len(BooP[((BooP.A==1)&(BooP.B==1))])
        X_10 = len(BooP[((BooP.A==1)&(BooP.B==0))])
        X_01 = len(BooP[((BooP.A==0)&(BooP.B==1))])
        X_00 = len(BooP[((BooP.A==0)&(BooP.B==0))])
        
        # Improvement (I) = pB/pA-1
        if X_11+X_10 == 0:
            I_x = 10101 # approx infinity!
        else:
            I_x = (X_11+X_01)/(X_11+X_10)-1
        
        # Results
        I_T[R] = I_x
    
    # Results
    I[BS] = I_T
    
    # CI
    Lower[BS] = np.percentile(list(I[BS].values()),  2.5)
    Media[BS] = np.percentile(list(I[BS].values()), 50  )
    Upper[BS] = np.percentile(list(I[BS].values()), 97.5)

print('Save')
output = open('MAE_B3_I.py.pkl', 'wb')
pickle.dump(I, output)
output.close()

output = open('MAE_B3_Lower.py.pkl', 'wb')
pickle.dump(Lower, output)
output.close()

output = open('MAE_B3_Media.py.pkl', 'wb')
pickle.dump(Media, output)
output.close()

output = open('MAE_B3_Upper.py.pkl', 'wb')
pickle.dump(Upper, output)
output.close()

print('Plot')
df_G1 = pd.DataFrame({'BS' : bootstrap_size,
                      'I'  : list(Media.values()),
                      'Lo' : list(Lower.values()),
                      'Hi' : list(Upper.values())})

fig, ax1 = plt.subplots(1,1)

df_G1.plot(x='BS', y='Hi', legend=False, ax=ax1, label='95% CI', linewidth=5, color='k', linestyle='--')
df_G1.plot(x='BS', y='I',  legend=False, ax=ax1, label='I',      linewidth=5, color='k', linestyle='-')
df_G1.plot(x='BS', y='Lo', legend=False, ax=ax1, label='95% CI', linewidth=5, color='k', linestyle='--')

for item in ([ax1.title, ax1.xaxis.label, ax1.yaxis.label] +
             ax1.get_xticklabels() + ax1.get_yticklabels()):
    item.set_fontsize(22)

legend = ax1.legend(loc=0, ncol=1, bbox_to_anchor=(0.9, -.3, 1, 1),
           fancybox=True, shadow=False,
           framealpha=1, fontsize=22) # , title='Percentile'
plt.setp(legend.get_title(),fontsize=22)

plt.xlabel('$BS$')
plt.ylabel('$I$')
plt.grid(b=True, which='major', color='b')
plt.grid(b=True, which='minor', color='b')
plt.xlim([0,BS_Max])
plt.ylim([0,2])

fig = plt.gcf()
fig.set_size_inches(4,4)
plt.show()
plt.clf()