# Probability in Python

# Importing packages

You can ignore this part for now.


```python
import numpy as np
```

# Import statistics module
We will use scipy.stats, which has several functions for statistics and probability distributions. 


```python
import scipy.stats as st
```

# Import pandas, matplotlib


```python
import pandas as pd
import matplotlib.pyplot as plt
```

# Function for uniform outcome

$n$: number of outcomes in the sample space

Output: $m$ outcomes selected uniformly at random from 1 to $n$


```python
def uniform(n, m):
  return np.random.randint(1, n+1, size = m)
```

# Toss a coin

Toss once, 10 times and 100 times

1: Heads and 2: Tails


```python
print(uniform(2, 1))
print(uniform(2, 10))
print(uniform(2,100))
```

    [1]
    [2 1 1 1 1 1 2 2 1 2]
    [1 1 2 2 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 1 1 2 2 2 1 1 1 2 2 2 1 2 1 1 2
     1 2 2 1 1 2 1 2 1 1 2 1 2 2 2 2 1 1 2 1 1 1 2 2 1 2 2 1 2 2 1 2 2 2 1 1 2
     1 2 1 1 2 1 1 2 1 2 1 1 2 1 2 2 1 2 1 2 1 1 1 1 1 2]


# Throw a die

Throw once, 10 times and 100 times


```python
print(uniform(6, 1))
print(uniform(6, 10))
print(uniform(6,100))
```

    [4]
    [6 2 2 5 3 2 3 1 1 3]
    [3 6 2 3 1 1 3 5 5 3 5 5 2 6 1 4 5 1 6 2 6 5 3 3 1 2 3 3 6 2 6 5 6 5 3 3 1
     2 3 3 1 6 1 5 6 5 5 2 4 3 4 6 2 5 6 2 2 3 3 5 2 1 1 4 3 2 1 1 4 6 1 4 2 6
     4 5 4 5 5 2 3 1 2 1 5 1 1 3 5 3 5 3 5 6 1 4 5 1 6 1]


# Estimating probability by simulation - Monte Carlo

The probability of an event $A$ can be estimated as follows. We can simulate the experiment repeatedly and independently, say $N$ times, and count the number of times the event occurred, say $N_A$. 

A good estimate of $P(A)$ is the following:
$$P(A) \approx \frac{N_A}{N}$$
As $N$ grows larger and larger, the estimate becomes better and better. This method is generally termed as Monte Carlo simulation.

We will first evaluate probability of coin toss described above using Monte Carlo simulations. There are two steps: generate a large number of tosses and count the number of heads or tails. These two steps can be written in a single loop usually.

You should run the simulation multiple times to see what probability estimate is obtained each time. You will see that the estimate is close to 0.5.


```python
no_heads = 0   #variable for storing number of heads
for i in range(1000): #repeat 1000 times
  if uniform(2, 1) == 1: #check if coin toss is heads
    no_heads = no_heads + 1
print(no_heads/1000) #probability estimate by Monte Carlo
```

    0.492


# Probability of die showing a number

We will modify the Monte Carlo simulation above for finding the probability that a dies shows a number falling in an event $A$. You will see that the estimate is close to $P(A)$. If you change the loop iterations to 10000, the estimate will be much closer to $P(A)$ and more consistent as well. 


```python
no = 0   #variable for storing number of event occurence
for i in range(10000): #repetitions
  die = uniform(6,1)  #experiment
  if die == 1 or die == 3: #Event
    no = no + 1
print(no/10000) #probability estimate by Monte Carlo
```

    0.3361


# Birthday problem

In a group of $n$ persons, what is the chance that some two have the same birthday? Assume birthday of a person is uniformly distributed in $\{1,2,\ldots,365\}$ and is independent of all other birthdays. Most people will think that you need at least 100 persons before you start seeing same birthdays. However, surprisingly perhaps, even with 23 persons there is a 50% chance of two sharing a birthday.

Event $A$: some two have same birthday

Event $A^c$: no two have same birthday

$A^c$: (Birthday 1 on any date $B_1$) and (Birthday 2 on any date other than $B_1$) and (Birthday 3 on any date other than $B_1$, $B_2$) and ... and (Birthday $n$ on any day other than $B_1,B_2,\ldots,B_{n-1}$)

$P(A^c)= 1 \cdot \left(1 - \frac{1}{365}\right)\left(1 - \frac{2}{365}\right)\cdots\left(1 - \frac{n-1}{365}\right)$

If $n=10$, what is the chance? If $n=30$, what is the chance?

We will do a Monte Carlo simulation to estimate the probability and compare with the calculation above.


```python
no = 0   #variable for storing number of event occurence
n = 60 #number of persons
print(1 - np.prod(1-np.arange(1,n)/365)) #probability from expression

for i in range(1000):
  B = np.zeros(366)   #array to keep track of birthdays seen
  for j in range(n): #generate birthdays for each person
    Bi = uniform(365, 1) #i-th birthday
    if B[Bi]  == 0:  #if Bi is seen for the first time
      B[Bi] = 1  #make note that Bi has been seen
    else:
      no = no + 1   #if Bi has been seen before, then two birthdays are same 
      break   #we can stop generating more birthdays and exit loop early

print(no/1000) #probability estimate by Monte Carlo
```

    0.994122660865348
    0.996


# Monty Hall problem

Here is the problem taken from the [Wiki page](https://en.wikipedia.org/wiki/Monty_Hall_problem).

> Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what's behind the doors, opens another door, say No. 3, which has a goat. He then says to you, "Do you want to pick door No. 2?" Is it to your advantage to switch your choice?

The assumptions (also taken from [Wiki](https://en.wikipedia.org/wiki/Monty_Hall_problem)) are as follows:
1. Car and goats are placed at random behind the doors.
2. Host always picks a door not chosen by contestant.
3. Host always reveals a goat and not a car.
4. Host always offers a choice to switch from the original door to the other closed door.

Under the above assumptions, here are the probabilities of winning.

P(win if contestant chooses to switch) = 2/3

P(win if contestant does not switch) = 1/3

You can see the Wiki page for the computation. Let us simulate and find the probability of winning under switch by Monte Carlo.


```python
no = 0   #variable for storing number of event occurence
for i in range(1000):
  car_loc = uniform(3, 1)
  if car_loc == 1:
    goat1_loc = 2; goat2_loc = 3
  elif car_loc == 2:
    goat1_loc = 1; goat2_loc = 3
  else:
    goat1_loc = 1; goat2_loc = 2

  contestant_orig = uniform(3, 1)
  if contestant_orig == goat1_loc:
    host_reveal_loc = goat2_loc; other_closed_door = car_loc
  elif contestant_orig == goat2_loc:
    host_reveal_loc = goat1_loc; other_closed_door = car_loc
  else:
    host_reveal_loc = goat1_loc; other_closed_door = goat2_loc
  if other_closed_door == car_loc:
    no = no + 1

print(no/1000) #probability estimate by Monte Carlo
```

    0.663


# Polya's urn scheme

Suppose an urn contains $r$ red and $b$ blue balls. The experiment proceeds in multiple steps, where Step $i$ is as follows:

Step $i$: Draw a ball at random, note down its colour and replace it in the urn. Add $c$ more balls of the same colour to the urn.

Let $R_i$ be the event that the $i$-th ball drawn is red. Let $B_i$ be the event that the $i$-th abll drawn is black.

Clearly, $P(R_1) = \frac{r}{r+b}$ and $P(B_1)=\frac{b}{r+b}$. It is perhaps surprising that, irrespective of $c$, we have, for all $i$,
$$P(R_i) = \frac{r}{r+b}, P(B_i) = \frac{b}{r+b}.$$
To prove the above, you can use induction. Assume that the above is true for $i$ and show it is true for $i+1$. Starting with $i=1$, by induction, the statement becomes true.

We will setup a Monte Carlo simulation for verifying $P(R_i)$ above for a few steps.


```python
no = 0   #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print(r/(r+b))
for i in range(1000):
  r = 10; b = 5
  c = 3
  for j in range(5): #do 5 steps
    if uniform(r+b, 1) <= r:
      r = r + c
    else:
      b = b + c
  if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
    no = no + 1
print(no/1000) #probability estimate by Monte Carlo
```

    0.6666666666666666
    0.66


# Gambler's ruin (simple random walk)

A gambler starting with $k$ units of money plays the following game at a casino:

* If he has $\ge 1$ units of money, a coin is tossed. If heads, the casino pays him 1 unit. If tails, he loses 1 unit to the casino.
* If he loses all money, he goes bankrupt and stops.
* If he gets $N$ units of money, he wins and stops playing.

If $p$ is the probability of heads and $q=1-p$, it can be shown that
$$\text{Pr}(\text{Bankruptcy})=\begin{cases}
1-k/N,&\text{ if }p=q=1/2,\\
\frac{\left(\dfrac{q}{p}\right)^k-\left(\dfrac{q}{p}\right)^N}{1-\left(\dfrac{q}{p}\right)^N}, &\text{ if }p\ne q.
\end{cases}$$
You can see some details of the proof of the above in the [Wiki page](https://en.wikipedia.org/wiki/Gambler%27s_ruin). Suppose $x_k$ denotes the probability of bankruptcy starting with $k$ units. The main idea is to condition on the first toss and derive the following recursive equation:
$$\begin{align}
x_k&=P(\text{Bankruptcy}\ |\ \text{first toss is head})\ p\ +\ P(\text{Bankruptcy}\ |\ \text{first toss is tail})\ q\\
&=x_{k+1}p+x_{k-1}q
\end{align}$$
with boundary conditions $x_0=1$ and $x_N=0$. Solution of the recursive equation results in the above closed form expression for $x_k$.

We are interested in Monte Carlo simulation of Gambler's ruin and verification of the formula for $x_k$. First, we consider the case $p=1/2$.


```python
no = 0   #variable for storing number of event occurence
k = 5; N = 10
print(1-k/N)
for i in range(1000):
  k = 5
  while k > 0 and k < N:
    if uniform(2, 1) == 1:
      k = k + 1
    else:
      k = k - 1
  if k == 0:
    no = no + 1
print(no/1000) #probability estimate by Monte Carlo
```

    0.5
    0.505


# Toss a biased coin

For $p\ne q$, we require a method to toss a biased coin. This is accomplished by the following function that generates $m$ coin tosses with probability of heads equal to $p$. Note that a value of 1 represents heads and 2 represents tails as before.


```python
def biased(p, m):
  return 2-(np.random.rand(m) < p)
```


```python
no_heads = 0   #variable for storing number of heads
p = 0.25
print(p)
for i in range(1000):
  if biased(p, 1) == 1:
    no_heads = no_heads + 1
print(no_heads/1000) #probability estimate by Monte Carlo
```

    0.25
    0.253


# Biased Gambler's ruin

We now simulate the biased version of Gambler's ruin.


```python
no = 0   #variable for storing number of event occurence
p = 0.35
qbyp = (1-p)/p
k = 5; N = 10
print((qbyp**k-qbyp**N)/(1-qbyp**N))
for i in range(1000):
  k = 5
  while k > 0 and k < N:
    if biased(p, 1) == 1:
      k = k + 1
    else:
      k = k - 1
  if k == 0:
    no = no + 1
print(no/1000) #probability estimate by Monte Carlo
```

    0.9566941509920124
    0.964


# Casino die game
Throw a pair of die. A player bets $k_1$ units of money on whether the sum of the two numbers is Under 7 or Over 7, and $k_2$ units on Equal to 7. For Under 7 and Over 7, the returns are $a$:1, while, for Equal to 7, the returns are $b$:1, if the player wins the bet. If the bet is lost, the unit of money goes to the casino.

The strategy for betting will be to independently and randomly select one of the 3 bets. The simulation will track the average return over a large number of trails.


```python
a = 1.0; b = 4.0
k1 = 1; k2 = 1
print((((a-1)*5-7)*k1+((b-1)-5)*k2)/6/3) #expected gain
avg_return = 0
for i in range(1000):
  bet = uniform(3,1) #1 - Under 7, 2 - Over 7, 3 - Equal to 7
  sum = uniform(6,1) + uniform(6,1)
  if ((bet == 1) and (sum < 7)) or ((bet == 2) and (sum > 7)): #win for Under 7 or Over 7 bet
    avg_return = avg_return + k1*(a-1)/1000
  if (bet == 3) and (sum == 7): #win for Equal to 7 bet
    avg_return = avg_return + k2*(b-1)/1000
  if ((bet == 1) and (sum >= 7)) or ((bet == 2) and (sum <= 7)): #loss for Under 7 or Over 7 bet
    avg_return = avg_return + (-k1)/1000
  if (bet == 3) and (sum != 7): #loss for Equal to 7 bet
    avg_return = avg_return + (-k2)/1000
  
print(avg_return) #simulated gain
```

    -0.5
    -0.5080000000000005


# Expected value of common distributions
The module has functions for generating binomial, geometric, Poisson and other distributions. We will generate a large number of samples and compute the average value and compare with the expected value.


```python
#binomial(20,0.3)
print(20*0.3) #expected value
x = st.binom.rvs(20,0.3,size=1000)
print(np.sum(x)/1000) #average value in simulation
```

    6.0
    6.062



```python
#geometric(0.3)
print(1/0.3) #expected value
x = st.geom.rvs(0.3,size=1000)
print(np.sum(x)/1000) #average value in simulation
```

    3.3333333333333335
    3.276



```python
#Poisson(6)
print(6) #expected value
x = st.poisson.rvs(6,size=1000)
print(np.sum(x)/1000) #average value in simulation
```

    6
    6.029


# Balls and bins
Suppose $m$ balls are thrown independently and uniformly at random into $n$ bins. We will compute the expected number of empty bins by simulation and compare with the theoretical value of $n(1-1/n)^m\approx ne^{-m/n}$.


```python
m = 10; n = 3
print(n*((1-1/n)**m)) #expected value
avg_empty_bins = 0
for i in range(1000):
  no_balls = np.zeros(n, dtype=int) #keep track of balls in bins
  for ball in range(m):
    bin = uniform(n, 1)
    no_balls[bin-1] += 1

  no_empty_bins = 0
  for bin in range(n):
    if no_balls[bin] == 0:
      no_empty_bins += 1

  avg_empty_bins += no_empty_bins/1000.0

print(avg_empty_bins) #average value in simulation
```

    0.0520245897474979
    0.05600000000000004


# Common continuous distributions and histograms
Scipy stats module can be used to generate samples from common continuous distributions. We will generate a number of such samples and plot their histogram to confirm that the samples follow the expected density function.

For histograms, we will use the hist() function from the matplotlib.pyplot module imported below.

## Uniform distribution
We will begin with the uniform distribution.


```python
x = st.uniform.rvs(0,3,size=10000)
plt.hist(x,bins=50,range=(0,3),density=True) #blue histogram
plt.plot([-0.2,0,0,3,3,3.2],[0,0,1.0/3,1.0/3,0,0],lw=2) #orange line. uniform[0,3] density
plt.show()
```

## From histogram to density
The code above generates 10000 samples that are supposed to be independent and uniformly distributed in $[0,3]$. The histogram, created using the plt.hist command, uses 100 bins of equal width in the range $[0,3]$. So, the bins are $[0,0.03),[0.03,0.06),\ldots,[2.97,3]$. 

Suppose the number of samples that fall into the bin $[0,0.03]$ is $N_0$. Then, by Monte Carlo, we have that 
$$P(0<X\le 0.03)\approx \frac{N_0}{10000},$$ 
where $X$ is a random variable with the sample distribution. Assuming that the density of $X$ satisfies $f_X(x)\approx f_X(0.015)$ over the bin ($0.015$ is the midpoint of the bin), we get
$$P(0<X\le 0.03)\approx 0.03f_X(0.015)\approx\frac{N_0}{10000}.$$
Using the above, we get 
$$f_X(0.015)\approx\frac{N_0}{300}.$$ 
Similarly, if $N_i$ is the number of samples in the $i$-th bin with midpoint $x_i$, we have
$$f_X(x_i)\approx\frac{N_i}{300}.$$
The option density=True in the plt.hist command specifies that the above calculation is to be done.

The plt.plot command plots the expected PDF as a line plot. The parameter lw specifies the linewidth and 2 pts.

Try changing the bin size to see how the plot changes. Does 50 bins look better? Why?

## Exponential distribution
We will next repeat the same for $X\sim$ Exp$(\lambda)$. The PDF is
$$f_X(x)=\lambda\exp(-\lambda x),$$
where $\lambda$ is called the scale parameter. Try changing the various parameters below to see what happens.


```python
x = st.expon.rvs(scale=1,size=10000)
plt.hist(x,bins=50,range=(0,10),density=True) #blue histogram
xp = np.linspace(0,10,50)
plt.plot(xp,st.expon.pdf(xp,scale=1),lw=2) #orange line, exp(\lambda) density
plt.show()
```

## Normal distribution
We will repeat the same for $X\sim$ Normal$(\mu,\sigma^2)$. The PDF is
$$f_X(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp(-(x-\mu)^2/2\sigma^2),$$
where the mean $\mu$ is called the location parameter `loc' and $\sigma$ is called the scale parameter.


```python
x = st.norm.rvs(loc=0, scale=1, size=10000)
plt.hist(x,bins=50,range=(-5,5),density=True) #blue histogram
xp = np.linspace(-5,5,50)
plt.plot(xp,st.norm.pdf(xp,loc=0,scale=1),lw=2) #orange line, normal pdf
plt.show()
```


```python
from sklearn.datasets import load_iris
```


```python
iris = load_iris()
print(iris.DESCR)
```

    .. _iris_dataset:
    
    Iris plants dataset
    --------------------
    
    **Data Set Characteristics:**
    
        :Number of Instances: 150 (50 in each of three classes)
        :Number of Attributes: 4 numeric, predictive attributes and the class
        :Attribute Information:
            - sepal length in cm
            - sepal width in cm
            - petal length in cm
            - petal width in cm
            - class:
                    - Iris-Setosa
                    - Iris-Versicolour
                    - Iris-Virginica
                    
        :Summary Statistics:
    
        ============== ==== ==== ======= ===== ====================
                        Min  Max   Mean    SD   Class Correlation
        ============== ==== ==== ======= ===== ====================
        sepal length:   4.3  7.9   5.84   0.83    0.7826
        sepal width:    2.0  4.4   3.05   0.43   -0.4194
        petal length:   1.0  6.9   3.76   1.76    0.9490  (high!)
        petal width:    0.1  2.5   1.20   0.76    0.9565  (high!)
        ============== ==== ==== ======= ===== ====================
    
        :Missing Attribute Values: None
        :Class Distribution: 33.3% for each of 3 classes.
        :Creator: R.A. Fisher
        :Donor: Michael Marshall (MARSHALL%PLU@io.arc.nasa.gov)
        :Date: July, 1988
    
    The famous Iris database, first used by Sir R.A. Fisher. The dataset is taken
    from Fisher's paper. Note that it's the same as in R, but not as in the UCI
    Machine Learning Repository, which has two wrong data points.
    
    This is perhaps the best known database to be found in the
    pattern recognition literature.  Fisher's paper is a classic in the field and
    is referenced frequently to this day.  (See Duda & Hart, for example.)  The
    data set contains 3 classes of 50 instances each, where each class refers to a
    type of iris plant.  One class is linearly separable from the other 2; the
    latter are NOT linearly separable from each other.
    
    .. topic:: References
    
       - Fisher, R.A. "The use of multiple measurements in taxonomic problems"
         Annual Eugenics, 7, Part II, 179-188 (1936); also in "Contributions to
         Mathematical Statistics" (John Wiley, NY, 1950).
       - Duda, R.O., & Hart, P.E. (1973) Pattern Classification and Scene Analysis.
         (Q327.D83) John Wiley & Sons.  ISBN 0-471-22361-1.  See page 218.
       - Dasarathy, B.V. (1980) "Nosing Around the Neighborhood: A New System
         Structure and Classification Rule for Recognition in Partially Exposed
         Environments".  IEEE Transactions on Pattern Analysis and Machine
         Intelligence, Vol. PAMI-2, No. 1, 67-71.
       - Gates, G.W. (1972) "The Reduced Nearest Neighbor Rule".  IEEE Transactions
         on Information Theory, May 1972, 431-433.
       - See also: 1988 MLC Proceedings, 54-64.  Cheeseman et al"s AUTOCLASS II
         conceptual clustering system finds 3 classes in the data.
       - Many, many more ...


# Summarizing the data
iris.data: 2D array containing all data (size 150 x 4)
 - Class 0: Row 0 to Row 49 of iris.data array
 - Class 1: Row 50 to Row 99 of iris.data array
 - Class 2: Row 100 to Row 149 of iris.data array

From the module scipy.stats, we can use the 'describe' command to get summary statistics of an array of data.


```python
# all data
sum_stats = st.describe(iris.data)
print(sum_stats.minmax)
```

    (array([4.3, 2. , 1. , 0.1]), array([7.9, 4.4, 6.9, 2.5]))



```python
print(sum_stats.mean)
print(sum_stats.variance)
```

    [5.84333333 3.05733333 3.758      1.19933333]
    [0.68569351 0.18997942 3.11627785 0.58100626]



```python
#Class 0
sum_stats = st.describe(iris.data[:50,:])
print(sum_stats.minmax)
print(sum_stats.mean)
print(sum_stats.variance)
```

    (array([4.3, 2.3, 1. , 0.1]), array([5.8, 4.4, 1.9, 0.6]))
    [5.006 3.428 1.462 0.246]
    [0.12424898 0.1436898  0.03015918 0.01110612]


# Plotting histograms to get a sense of data
For every class, we can create histogram plots. The code below generates these histograms for Class 0. 




```python
plt.subplot(221)
plt.hist(iris.data[:50,0])
plt.xlim([0,6])
plt.title('SL')
plt.subplot(222)
plt.hist(iris.data[:50,1])
plt.xlim([0,6])
plt.title('SW')
plt.subplot(223)
plt.hist(iris.data[:50,2])
plt.xlim([0,6])
plt.title('PL')
plt.subplot(224)
plt.hist(iris.data[:50,3])
plt.xlim([0,6])
plt.title('PW')
plt.suptitle('Class 0')
plt.tight_layout()
plt.show()
```

# 2D histograms
The following code shows how to plot 2D histograms. This is being provided as a reference for those who are interested.


```python
fig = plt.figure(figsize=(10,5), dpi = 150)
ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122, projection='3d')

H, _x, _y = np.histogram2d(iris.data[:50,0], iris.data[:50,1])
_xx, _yy = np.meshgrid(_x[1:], _y[1:])
xx, yy, HH = _xx.ravel(), _yy.ravel(), H.ravel()

ax1.bar3d(xx, yy, np.zeros_like(HH), 0.13, 0.13, HH, shade=True, color='b')

ax1.set_xlabel('SL')
ax1.set_ylabel('SW')

H, _x, _y = np.histogram2d(iris.data[:50,2], iris.data[:50,3])
_xx, _yy = np.meshgrid(_x[1:], _y[1:])
xx, yy, HH = _xx.ravel(), _yy.ravel(), H.ravel()

ax2.bar3d(xx, yy, np.zeros_like(HH), 0.04, 0.04, HH, shade=True, color='g')

ax2.set_xlabel('PL')
ax2.set_ylabel('PW')
```

# Fitting a distribution and estimating parameters
Interarrival times in a Gamma ray experiment are given in the book "Mathematical Statistics and Data Analysis" by John Rice. The data is in the form of an excel file. 

We will use the pandas module in python to read the excel file.


```python
#Make sure you upload the file gamma-arrivals.xls to the Colab runtime before running this...
df = pd.read_excel(r'/content/gamma-arrivals.xls')
df
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Col1</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>38.026635</td>
    </tr>
    <tr>
      <th>1</th>
      <td>409.445210</td>
    </tr>
    <tr>
      <th>2</th>
      <td>272.463390</td>
    </tr>
    <tr>
      <th>3</th>
      <td>98.674526</td>
    </tr>
    <tr>
      <th>4</th>
      <td>87.858008</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
    </tr>
    <tr>
      <th>3930</th>
      <td>321.296880</td>
    </tr>
    <tr>
      <th>3931</th>
      <td>146.966750</td>
    </tr>
    <tr>
      <th>3932</th>
      <td>184.342160</td>
    </tr>
    <tr>
      <th>3933</th>
      <td>84.361201</td>
    </tr>
    <tr>
      <th>3934</th>
      <td>218.148090</td>
    </tr>
  </tbody>
</table>
<p>3935 rows × 1 columns</p>
</div>



We see that there are 3935 samples of data. To get a sense of the distribution, we should plot a histogram.


```python
plt.hist(df['Col1'],bins=50)
```

# Fitting a Gamma distribution
From the histogram, the distribution could be modelled as Gamma$(\alpha,\beta)$. The next step is to estimate $\alpha$ and $\beta$ from the given samples.

## Method of moments
Suppose $m_1$ and $m_2$ are the first and second moments of the samples. The method of moments estimates are obtained by solving
$$m_1=\frac{\alpha}{\beta},$$
$$m_2=\frac{\alpha^2}{\beta^2}+\frac{\alpha}{\beta^2}.$$
The solution results in
$$\hat{\alpha}_{MM}=\frac{m_1^2}{m_2-m_1^2}=\frac{m_1^2}{s^2},\hat{\beta}_{MM}=\frac{m_1}{m_2-m_1^2}=\frac{m_1}{s^2}.$$
We now compute the values of $m_1$ (sample mean) and $s^2=m_2-m_1^2$ (sample variance) from the data. After that, we can compute the estimates.


```python
x = np.array(df['Col1'])
m1 = np.average(x)
ss = np.var(x)
print(m1)
print(ss)

```

    79.93522014279034
    6311.67624499251



```python
alphaMM = m1*m1/ss
betaMM = m1/ss
print(alphaMM)
print(betaMM)
```

    1.0123522137792949
    0.012664657856335469


We can plot the density of the Gamma on top of the density histogram to check if the estimate gives a reasonable fit.


```python
fig,ax = plt.subplots(1,1)
ax.hist(x,density=True,bins=50)
xx = np.linspace(0,300,50)
ax.plot(xx, st.gamma.pdf(xx,alphaMM,scale=1/betaMM),label='gamma fit MM')
ax.legend(loc='best')
plt.show()
```

# Bootstrap
How do we find the bias and variance of the estimator? Theoretical derivations of the sampling distributions may be too cumbersome and difficult in most cases. Bootstrap is a Monte Carlo simulation method for computing metrics such as bias, variance and confidence intervals for estimators.

In the above example, we have found $\hat{\alpha}_{MM}=1.0123...$ and $\hat{\beta}_{MM}=0.01266...$. Using these values, we simulate $n=3935$ *iid* samples from Gamma$(1.0123...,0.0126...)$ and, using the simulated samples, we compute new estimates of $\alpha$ and $\beta$ and call them $\hat{\alpha}_{MM}(1)$ and $\hat{\beta}_{MM}(1)$. Now, repeat the simulation $N$ times to get estimates $\hat{\alpha}_{MM}(i)$ and $\hat{\beta}_{MM}(i)$, $i=1,2,\ldots,N$.

The sample variance of $\{\hat{\alpha}_{MM}(1), \hat{\alpha}_{MM}(2),\ldots,\hat{\alpha}_{MM}(N)\}$ is taken to be the bootstrap estimate for the variance of the estimator.


```python
N = 1000
n = 3935
alpha_hat = np.zeros(N)
beta_hat = np.zeros(N)
for i in np.arange(N):
  xi = st.gamma.rvs(alphaMM,scale=1/betaMM,size=n)
  m1i = np.average(xi); ssi = np.var(xi)
  alpha_hat[i] = m1i*m1i/ssi; beta_hat[i] = m1i/ssi
```

We can see the histograms of the estimates to get an idea of the spread of the values.


```python
ax1 = plt.subplot(121)
ax1.hist(alpha_hat,density=True)
ax2 = plt.subplot(122)
ax2.hist(beta_hat,density=True)
```

Notice how the histograms look roughly normal.

The sample standard deviations of the estimates is a bootstrap estimate for the standard error of the estimator.


```python
print(np.sqrt(np.var(alpha_hat)))
print(np.sqrt(np.var(beta_hat)))
```

    0.031822357311654136
    0.000449581400441715


## Confidence intervals
Suppose a parameter $\theta$ is estimated as $\hat{\theta}$, and suppose the distribution of $\hat{\theta}-\theta$ is known. Then, to obtain $(100(1-\alpha))$% confidence intervals (typical values are $\alpha=0.1$ for 90% confidence intervals and $\alpha=0.05$ for 95% confidence intervals), we use the CDF of $\hat{\theta}-\theta$ to obtain $\delta_1$ and $\delta_2$ such that
$$P(\hat{\theta}-\theta\le\delta_1)=1-\frac{\alpha}{2},$$
$$P(\hat{\theta}-\theta\le\delta_2)=\frac{\alpha}{2}.$$
Actually, the inverse of the CDF of $\hat{\theta}-\theta$ is used to find the above $\delta_1$ and $\delta_2$. From the above, we see that
$$P(\hat{\theta}-\theta \le \delta_1)-P(\hat{\theta}-\theta \le \delta_2)= P(\delta_2< \hat{\theta}-\theta \le \delta_1)=1-\frac{\alpha}{2}-\frac{\alpha}{2}=1-\alpha.$$
The above is rewritten as
$$P(\hat{\theta}-\delta_1\le\theta<\hat{\theta}-\delta_2)=1-\alpha,$$
and $[\hat{\theta}-\delta_1,\hat{\theta}-\delta_2]$ is interpreted as the $100(1-\alpha)$% confidence interval.

## Bootstrap confidence intervals
The CDF of $\hat{\theta}-\theta$ might be difficult to determine in many cases, and the bootstrap method is used often to estimate $\delta_1$ and $\delta_2$. We consider the list of numbers $\{\hat{\alpha}_{MM}(1)-1.0123...,\ldots,\hat{\alpha}_{MM}(N)-1.0123...\}$ and pick the $100(\alpha/2)$-th percentile and $100(1-\alpha/2)$-th percentile.


```python
del1 = np.percentile(alpha_hat - alphaMM, 97.5)
del2 = np.percentile(alpha_hat - alphaMM, 2.5)
print([del1,del2])
```

    [0.06149537903410332, -0.06043329601919404]


The 95% confidence interval for $\alpha$ using the method of moments estimator works out to $[1.0123-0.0615,1.0123-(-0.0604)]=[0.9508,1.0727]$.

## Maximum likelihood
We now turn to the maximum likelihood estimator for $\alpha$ and $\beta$. The likelihood $L(x_1,\ldots,x_n)$ can be written as
$$L = \frac{\beta^\alpha}{\Gamma(\alpha)}x_1^{\alpha-1}e^{-\beta x_1}\,\frac{\beta^\alpha}{\Gamma(\alpha)}x_2^{\alpha-1}e^{-\beta x_2}\cdots \frac{\beta^\alpha}{\Gamma(\alpha)}x_n^{\alpha-1}e^{-\beta x_n}= \frac{\beta^{n\alpha}}{\Gamma(\alpha)^n}(x_1\cdots x_n)^{\alpha-1}e^{-\beta(x_1+\cdots+x_n)},$$
$$\log L = n\alpha\log\beta-n\log\Gamma(\alpha)+(\alpha-1)\log(x_1\cdots x_n)-\beta(x_1+\cdots+x_n).$$
Differentiating $\log L$ with respect to $\beta$ and equating to zero, we get
$$n\alpha\frac{1}{\beta}-(x_1+\cdots+x_n)=0,\text{or }\alpha=\beta \frac{x_1+\cdots+x_n}{n}.$$
Differentiating $\log L$ with respect to $\alpha$ and equating to zero, we get
$$n\log\beta-n\frac{\Gamma'(\alpha)}{\Gamma(\alpha)}+\log(x_1\cdots x_n)=0.$$
So, we get two equations in the two variables $\alpha$ and $\beta$. However, the equations do not have a closed form solution, and we need to solve them numerically or approximately. From the first equation, we have $\log\beta=\log\alpha-\log\frac{x_1+\cdots+x_n}{n}$. Using this in the second equation, we get
$$\log\alpha - \frac{\Gamma'(\alpha)}{\Gamma(\alpha)}=\log\frac{x_1+\cdots+x_n}{n}-\frac{1}{n}\log(x_1\cdots x_n).$$
We will now solve the above equation to find the ML estimate of $\alpha$. This will be a numerical solution.


```python
lm1 = np.average(np.log(x))
#Write the equation as a function
#digamma is the function Gamma'/Gamma
from scipy.special import digamma
fML = lambda a: (np.log(a) - digamma(a) - np.log(m1)+lm1)
```

We can plot the above function to see how it looks.


```python
fig, ax = plt.subplots(1,1)

xx = np.linspace(0.1,2,50)
ax.plot(xx,fML(xx))
ax.grid(True)
plt.show()
```


```python
#For solving numerically, we will use scipy.optimize
import scipy.optimize as sopt
sol = sopt.root_scalar(fML, bracket=[0.1,2])
sol.root
```




    1.0263317358725452




```python
alphaML = sol.root
betaML = alphaML/m1
print([alphaML, betaML])
```

    [1.0263317358725452, 0.012839543495835532]


Let us check the fit with the histogram.


```python
fig,ax = plt.subplots(1,1)
ax.hist(x,density=True,bins=50)
xx = np.linspace(0,300,50)
ax.plot(xx, st.gamma.pdf(xx,alphaMM,scale=1/betaMM),lw='4',alpha=0.7,label='gamma fit MM')
ax.plot(xx, st.gamma.pdf(xx,alphaML,scale=1/betaML),lw='1',label='gamma fit ML')
ax.legend(loc='best')
plt.show()
```

Both the curves are literally on top of each other showing very good fit. Let us use the bootstrap method to find variance and confidence intervals for the ML estimator.


```python
N = 1000
n = 3935
alpha_hatML = np.zeros(N)
beta_hatML = np.zeros(N)
for i in np.arange(N):
  xi = st.gamma.rvs(alphaMM,scale=1/betaMM,size=n)
  m1i = np.average(xi); lm1i = np.average(np.log(xi))
  fMLi = lambda a: (np.log(a) - digamma(a) - np.log(m1i)+lm1i)
  soli = sopt.root_scalar(fMLi, bracket = [0.1,2]) 
  alpha_hatML[i] = soli.root; beta_hatML[i] = soli.root / m1i
```


```python
print(np.sqrt(np.var(alpha_hatML)))
print(np.sqrt(np.var(beta_hatML)))
```

    0.01988186514427374
    0.0003157550854393247


We see that the variance of the bootstrap ML estimator is lesser than that of bootstrap MM estimator.

8188dcb4d1a510edd1b40f2d0a9de3be87020a28
ipynb
Jupyter Notebook
Statistics/stats.ipynb
tanav2202/100-days-python-and-Ml
70ac964a98103bc68ebdf1d03bebc625ace93cac
MIT
2021-12-17T19:47:00.000Z
2022-01-13T16:19:54.000Z
Statistics/stats.ipynb
tanav2202/100-days-python-and-Ml
70ac964a98103bc68ebdf1d03bebc625ace93cac
MIT
Statistics/stats.ipynb
tanav2202/100-days-python-and-Ml
70ac964a98103bc68ebdf1d03bebc625ace93cac
MIT
Qwen/Qwen-72B
1. YES
2. YES