<jupyter_start><jupyter_text><!--  TODO: Do bootstrap regression example fromStatistical_Models_Freedman.pdf -->As we have seen, it can be very difficultor impossible to determinethe probability density distribution of the estimatorof somequantity. The idea behind the bootstrap is that we can use computationto approximate  these functions which would otherwise be impossible tosolve foranalytically. Let's start with a simple example. Suppose we have the followingset of randomvariables, $\lbrace X_1, X_2, \ldots, X_n \rbrace$ where each $X_k\sim F$. Inother words the samples are all drawn from the same unknowndistribution $F$.Having run the experiment, we thereby obtain the followingsample set:$$\lbrace x_1, x_2, \ldots, x_n \rbrace$$ The sample mean is computed from this set as,$$\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$$ The next question is how close is the sample mean to the true mean,$\theta =\mathbb{E}_F(X)$. Note that the second central moment of $X$ is asfollows:$$\mu_2(F) := \mathbb{E}_F (X^2)  - (\mathbb{E}_F (X))^2$$ The standard deviation of the sample mean, $\bar{x}$, given $n$samples from anunderlying distribution $F$, is the following:$$\sigma(F) = (\mu_2(F)/n)^{1/2}$$ Unfortunately, because we have only the set of samples $\lbrace x_1,x_2,\ldots, x_n \rbrace$ and not $F$ itself, we cannot compute this andinstead mustuse the estimated standard error,$$\bar{\sigma} = (\bar{\mu}_2/n)^{1/2}$$ where $\bar{\mu}_2 = \sum (x_i -\bar{x})^2/(n-1)$, which is theunbiasedestimate of $\mu_2(F)$. However, this is not the only way to proceed.Instead,we could replace $F$ by some estimate, $\hat{F}$ obtained as apiecewisefunction of $\lbrace x_1, x_2, \ldots, x_n \rbrace$ by placingprobability mass$1/n$ on each $x_i$. With that in place, we can compute theestimated standarderror as the following:$$\hat{\sigma}_B  = (\mu_2(\hat{F})/n)^{1/2}$$ which is called the *bootstrap estimate* of the standard error.Unfortunately,the story effectively ends here. In even a slightly more generalsetting, thereis no clean formula $\sigma(F)$ within which $F$ can be swappedfor $\hat{F}$.This is where the computer saves the day. We actually do not need to know theformula $\sigma(F)$ because we can compute it using a resampling method. Thekeyidea is to sample with replacement from $\lbrace x_1, x_2, \ldots, x_n\rbrace$.The new set of $n$ independent draws (with replacement) from this setis the*bootstrap sample*,$$y^* = \lbrace x_1^*, x_2^*, \ldots, x_n^* \rbrace$$The Monte Carlo algorithm proceeds by first by selecting a large number  ofbootstrap samples, $\lbrace y^*_k\rbrace$, then computing the statistic on eachof these samples, and then computing the sample standard deviation of theresults in the usual way. Thus, the bootstrap estimate of the statistic$\theta$is the following,$$\hat{\theta}^*_B = \frac{1}{B} \sum_k \hat{\theta}^*(k)$$ with the corresponding square of the sample standard deviation as$$\hat{\sigma}_B^2 = \frac{1}{B-1} \sum_k (\hat{\theta}^*(k)-\hat{\theta}^*_B)^2$$ The process is much simpler than  the notation implies.Let's explore this witha simple example using Python. The next blockof code sets up some samples froma $\beta(3,2)$ distribution,<jupyter_code>import numpy as np
_=np.random.seed(123456)
import numpy as np
from scipy import stats
rv = stats.beta(3,2)
xsamples = rv.rvs(50)<jupyter_output><empty_output><jupyter_text>Because this is simulation data, we already know that themean is $\mu_1 = 3/5$and the standard deviation of the sample meanfor $n=50$ is $\bar{\sigma}=\sqrt{2}/50$, which we will verifylater.<jupyter_code>%matplotlib inline

from matplotlib.pylab import subplots
fig,ax = subplots()
fig.set_size_inches(8,4)
_=ax.hist(xsamples,density=True,color='gray')
ax2 = ax.twinx()
_=ax2.plot(np.linspace(0,1,100),rv.pdf(np.linspace(0,1,100)),lw=3,color='k')
_=ax.set_xlabel('$x$',fontsize=28)
_=ax2.set_ylabel('  $y$',fontsize=28,rotation='horizontal')
fig.tight_layout()
fig.savefig('fig-statistics/Bootstrap_001.png')<jupyter_output><empty_output><jupyter_text><!-- dom:FIGURE: [fig-statistics/Bootstrap_001.png, width=500 frac=0.85] The$\beta(3,2)$ distribution and the histogram that approximates it. <divid="fig:Bootstrap_001"> --><divid="fig:Bootstrap_001">The $\beta(3,2)$ distribution and thehistogram that approximates it.<img src="fig-statistics/Bootstrap_001.png"width=300>[Figure](fig:Bootstrap_001) shows the$\beta(3,2)$ distribution andthe corresponding histogram of the samples. Thehistogram represents$\hat{F}$ and is the distribution we sample from to obtainthebootstrap samples. As shown, the $\hat{F}$ is a pretty crude estimateforthe $F$ density (smooth solid line), but that's not a seriousproblem insofar asthe following bootstrap estimates are concerned.In fact, the approximation$\hat{F}$ has a naturally tendency topull towards the bulk of probability mass.This is afeature, not a bug; and is the underlying mechanism that explainsbootstrapping, but the formal proofs that exploit this basicidea are far out ofour scope here.  The next block generates thebootstrap samples<jupyter_code>yboot = np.random.choice(xsamples,(100,50))
yboot_mn = yboot.mean()<jupyter_output><empty_output><jupyter_text>and the bootstrap estimate is therefore,<jupyter_code>np.std(yboot.mean(axis=1)) # approx sqrt(1/1250)<jupyter_output><empty_output><jupyter_text>[Figure](fig:Bootstrap_002) shows the distribution of computedsample meansfrom the bootstrap samples.  As promised, the next blockshows how to use`sympy.stats` to compute the $\beta(3,2)$ parameters we quotedearlier.<jupyter_code>fig,ax = subplots()
fig.set_size_inches(8,4)
_=ax.hist(yboot.mean(axis=1),density=True,color='gray',rwidth=.8)
_=ax.set_title('Bootstrap std of sample mean %3.3f vs actual %3.3f'%
             (np.std(yboot.mean(axis=1)),np.sqrt(1/1250.)))
fig.tight_layout()
fig.savefig('fig-statistics/Bootstrap_002.png')<jupyter_output><empty_output><jupyter_text><!-- dom:FIGURE: [fig-statistics/Bootstrap_002.png, width=500 frac=0.85] Foreach bootstrap draw, we compute the sample mean. This is the histogram of thosesample means that will be used to compute the bootstrap estimate of the standarddeviation.  --><divid="fig:Bootstrap_002">For each bootstrap draw, we compute the samplemean. This is the histogram of those sample means that will be used to computethe bootstrap estimate of the standard deviation.<img src="fig-statistics/Bootstrap_002.png" width=300><jupyter_code>import sympy as S
import sympy.stats
for i in range(50): # 50 samples
    # load sympy.stats Beta random variables
    # into global namespace using exec
    execstring = "x%d = S.stats.Beta('x'+str(%d),3,2)"%(i,i)
    exec(execstring) 

# populate xlist with the sympy.stats random variables
# from above
xlist = [eval('x%d'%(i)) for i in range(50) ]
# compute sample mean
sample_mean = sum(xlist)/len(xlist)
# compute expectation of sample mean
sample_mean_1 = S.stats.E(sample_mean).evalf()
# compute 2nd moment of sample mean
sample_mean_2 = S.stats.E(S.expand(sample_mean**2)).evalf()
# standard deviation of sample mean
# use sympy sqrt function
sigma_smn = S.sqrt(sample_mean_2-sample_mean_1**2) # sqrt(2)/50
print(sigma_smn)<jupyter_output>0.0282842712474623<jupyter_text>**Programming Tip.**Using the `exec` function enables the creation of asequence of Sympyrandom variables. Sympy has the `var` function which canautomaticallycreate a sequence of Sympy symbols, but there is no correspondingfunction in the statistics module to do this for random variables.<!--@@@CODE src-statistics/Bootstrap.py from-to:^import sympy as S@^print(sigma_smn) -->**Example.** Recall the delta methodfrom the section [sec:delta_method](sec:delta_method).  Suppose we have a setof Bernoulli coin-flips($X_i$) with probability of head $p$. Our maximumlikelihood estimatorof $p$ is $\hat{p}=\sum X_i/n$ for $n$ flips.  We know thisestimatoris unbiased with $\mathbb{E}(\hat{p})=p$ and $\mathbb{V}(\hat{p}) =p(1-p)/n$. Suppose we want to use the data to estimate the variance oftheBernoulli trials ($\mathbb{V}(X)=p(1-p)$). By the notation thedelta method,$g(x) = x(1-x)$. By the plug-in principle, our maximumlikelihood estimator ofthis variance is then $\hat{p}(1-\hat{p})$. Wewant the variance of thisquantity.  Using the results of the deltamethod, we have$$\begin{align*}\mathbb{V}(g(\hat{p})) &=(1-2\hat{p})^2\mathbb{V}(\hat{p})\\\\mathbb{V}(g(\hat{p})) &=(1-2\hat{p})^2\frac{\hat{p}(1-\hat{p})}{n} \\\\end{align*}$$ Let's see how useful this is with a short simulation.<jupyter_code>import numpy as np
np.random.seed(123)
from scipy import stats
import numpy as np
p= 0.25 # true head-up probability
x = stats.bernoulli(p).rvs(10)
print(x)<jupyter_output>[0 0 0 0 0 0 1 0 0 0]<jupyter_text>The maximum likelihood estimator of $p$ is $\hat{p}=\sum X_i/n$,<jupyter_code>phat = x.mean()
print(phat)<jupyter_output>0.1<jupyter_text>Then, plugging this into the delta method approximant above,<jupyter_code>print((1-2*phat)**2*(phat)**2/10)<jupyter_output>0.0006400000000000003<jupyter_text>Now, let's try this using the bootstrap estimate of the variance<jupyter_code>phat_b=np.random.choice(x,(50,10)).mean(1)
print(np.var(phat_b*(1-phat_b)))<jupyter_output>0.0050490000000000005<jupyter_text>This shows that the delta method's estimated varianceis different from thebootstrap method, but which one is better?For this situation we can solve forthis directly using Sympy<jupyter_code>import sympy as S
from sympy.stats import E, Bernoulli
xdata =[Bernoulli(i,p) for i in S.symbols('x:10')]
ph = sum(xdata)/float(len(xdata))
g = ph*(1-ph)<jupyter_output><empty_output><jupyter_text>**Programming Tip.**The argument in the `S.symbols('x:10')` function returns asequence of Sympysymbols named `x1,x2` and so on. This is shorthand forcreating and naming eachsymbol sequentially. Note that `g` is the$g(\hat{p})=\hat{p}(1- \hat{p})$ whose variance we are trying to estimate.Then,we can plug in for the estimated $\hat{p}$ and get the correctvalue forthe variance,<jupyter_code>print(E(g**2) - E(g)**2)<jupyter_output>0.00442968750000000<jupyter_text>This case is generally representative --- the delta method tendstounderestimate the variance and the bootstrap estimate is better here. Parametric BootstrapIn the previous example, we used the $\lbrace x_1, x_2,\ldots, x_n \rbrace $samples themselves as the basis for $\hat{F}$ by weightingeach with $1/n$. Analternative is to *assume* that the samples come from aparticulardistribution, estimate the parameters of that distribution from thesample set,and then use the bootstrap mechanism to draw samples from theassumeddistribution, using the so-derived parameters. For example, the nextcode blockdoes this for a normal distribution.<jupyter_code>rv = stats.norm(0,2)
xsamples = rv.rvs(45)
# estimate mean and var from xsamples
mn_ = np.mean(xsamples)
std_ = np.std(xsamples)
# bootstrap from assumed normal distribution with
# mn_,std_ as parameters
rvb = stats.norm(mn_,std_) #plug-in distribution
yboot = rvb.rvs(1000)<jupyter_output><empty_output><jupyter_text>Recallthe sample variance estimator is the following:$$S^2 = \frac{1}{n-1} \sum (X_i-\bar{X})^2$$ Assuming that the samples are normally distributed, thismeans that$(n-1)S^2/\sigma^2$ has a chi-squared distribution with$n-1$ degrees offreedom. Thus, the variance, $\mathbb{V}(S^2) = 2\sigma^4/(n-1) $. Likewise,the MLE plug-in estimate for this is$\mathbb{V}(S^2) = 2 \hat{\sigma}^4/(n-1)$The following code computesthe variance of the sample variance, $S^2$ using theMLE and bootstrapmethods.<jupyter_code># MLE-Plugin Variance of the sample mean 
print(2*(std_**2)**2/9.)         # MLE plugin
# Bootstrap variance of the sample mean 
print(yboot.var())
# True variance of sample mean 
print(2*(2**2)**2/9.)<jupyter_output>2.22670148617726
3.2946788568183387
3.5555555555555554
{'fraction_non_alphanumeric': 0.1258989762247229, 'fraction_numerical': 0.033166934596835604, 'mean_word_length': 5.848203939745075, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '331754', 'license': 'MIT', 'n_tokens_mistral': 4066, 'n_tokens_neox': 3761, 'n_words': 1251, 'path': 'chapter/statistics/Bootstrap.ipynb', 'repo_name': 'mhdella/Python-for-Probability-Statistics-and-Machine-Learning-2E', 'size': '11819'}