---
title: "Combinatorics of the bootstrap"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---


# Motivation

   + bootstrapping combinatorial properties?

   + When we draw n items with replacement from a sample of size n, we may repeat samples.
   
    As a consequence, any bootstrap sample that contains repeats ignore a part of the original data.
 

# Elementary maths for the bootstrap

How do we calculate approximately $log(1+x)$ when $x$ is small?


Applying the Taylor formula near 1 gives:

$$\log(1+x)=x + \ldots $$

OK. Let's use R to check the quality of this approximation:

```{r}
n=30
for(n in 1:100){
x=1/n
print(log(1+x) - x)
}


```

How do we calculate the exponential of a number?

$$e^x = \lim_{n \rightarrow \infty} (1+\frac{x}{n})^{-n}$$


$$x = -n \log (1 + \frac{x}{n} )  =    x \frac{x}{n} + \ldots $$

OK. Let's use R to check the quality of this approximation:


```{r}
n=30
for(n in 1:100){
print(exp(1) - (1-1/n)^(-n))
}
```

To calculate the number of permutations of n items, we need to calculate the factorial

But this becomes very complicated for large values of $n$. Is there a quicker and simpler way to get an approximate value of factorial $n!$?

Yes there is: Stirling's formula. In , Stirling showed that:

$$n! \approx n^n e^{-n} \sqrt{2\pi n}$$



OK. Let's use R to check the quality of this approximation:

```{r}
e=exp(1)
for(n in 1:20){
sn=(n^n) * e^(-n) * (2 *pi * n)^(0.5)
print(c(factorial(n),sn,(factorial(n) - sn)/sn))
}
```

How would you use S to calculate approximate values of the binomial coefficients?


Remember that the number of different sets of $p$ items chosen among $n$ is:

$$\binom{m}{n} = \frac{n!}{p!(n-p)!}$$


# How many samples does bootstrap samples ignore on average?



Suppose you roll six identical six-sided dice. 
Chance are that you will see at least one repeated number.
 
The probability that you will see six unique numbers is 
 
   
 ```{r eval=F}
set.seed(117426)

myListOfStats = NULL
nx=6
B=10
ip=1
is = 0
cval <- c(0,0,0,0,0,0)

data = 1:nx
for (i in 1:B){
xb   =  sample(data,nx,replace=T)
#print(table(xb))
for(k in 1:nx){
    if(is == 1){
	cval[k] = cval[k] + sum(xb == k)
	}
	if(is == 0){
	cval[k] =  sum(xb == k)
	}
}
if(ip == 1){
print(xb)
print(cval)
print('------------')
}
}

#mean(cval/B)
```
 
   
Recall the  "birthday matching problem" : what is the probability that 2 persons (at least) in a group of n persons 
have their birthday on the same day?
   
This is more than 50% for $n>=23$.
  
avoid u and avoid u and avoid u


$$p=(1-1/n)^n$$


Apply our maths. When n is large:

$$p = 1/e \approx 0.368$$


That is how much of the original data will be ignored on average in bootstrap samples.


```{r}
set.seed(117426)
nx=100 # try 8 15 30 100
B=100000
# then try B=1000000 no diff for nx=100
B=1000000
cval <- matrix(0,ncol=nx,nrow=1)

data = 1:nx
for (i in 1:B){
xb   =  sample(data,nx,replace=T)
#print(table(xb))
for(k in 1:nx){
if(k %in% xb){cval[k]=cval[k]+1}
   # if(sum(xb == k) == 0) cval[k] = cval[k]+1
}
#print(cval)
}

mean(cval/B) + exp(-1)

```

On average the sample contains 63.2% of the original observations and omits 36.8%. 
   
   In other words, the average bootstrap sample omits 36.8% of the original data.
   
   What about the variability of this proportion?
   
   This fraction $p_m$ will be approximately normally distributed in the large sample limit 
   
   $$p_m \sim N(1/e, \sigma^2 = p(1-p)/n)$$


## is it a problem that many bsamples do not see a big chunck of the original data?


That is how much of the original data will be ignored on average in bootstrap samples.

Is this a problem


# Counting the unique patterns produced by the bootstrap


Exchangeable statistics (sample min: Notexchangeable, sample mean: exchangeable)

Do not distinguish two sequences that are permutations of each other


Let us try to count unique patterns

METHOD 1

Recoding to avoid redundancy: let us call 'pattern' the frequency of original data in the bootstrap sample:

11111    [50000]

32233    [02300]     33322  32332     


hence in the code space we count unique patterns !

coding = a vector of counts where each count can take values between 0 and n

constraint is : 

$$\sum_{i=1}^n m_i = n$$


METHOD 2

Allocate n symbols to n boxes. THERE MUST BE n BOXES EVEN if SOME BOXES STAY EMPTY

the BOXES are the ORIGINAL LABELS of the data

How many patterns should we expect. Let us again find another way to code the patterns.

I will represent a pattern by a symbol 

$$212113 =  111 | 22 | 3| |$$


I can also use the coding:

How many ways to we have to place $n-1$ separators in $n+n-1$ places.

For 4 objects ($n=4$)


$$C_n = \binom{2n-1}{n-1}$$

$$C_n \approx (n \pi)^{-1/2}2^{2n-1}$$

$$n! \approx n^n e^{-n} \sqrt{2\pi n}$$

AND NOW we can apply our exercises. Do this yourself, calculate $C_n$ approximately


```{r}
# your R code here

```

# Which pattern is going to be the most represented?

Can you test your presumption with an R code?


```{r}
# your R code here

```


# Approximate normality of boostrap distributions



```{r}
set.seed(117426)
dof=0.2
n=200000 # 20 200 2000 20000 200000
data=rchisq(n,dof,ncp=5)
median(data)
mean(data)
B=10000 # 100000 1000000
bmedian = NULL
for (i in 1:B){
mybsample   =  sample(data,n,replace=T)
#mystat  =  mystatchoice(1, mybsample)
mystat = median(mybsample)
bmedian = c(bmedian,mystat)
}
hist(bmedian)
```


What should we expect?

$$\sigma^2 = \frac{1}{4 f(m)^2}$$


Links to sources and online documents:

<ul>
<li><a href="https://blogs.sas.com/content/iml/2017/06/28/average-bootstrap-sample-omits-data.htm">Blog : average-bootstrap-sample-omits-data</a>
<li><a href="http://statweb.stanford.edu/~susan/courses/s208/node37.html">Susan Holmes's lectures at Stanford</a>
</ul>
