#### Param Estimation

* We randomly count the number of observed moving traffic citations within a 5 mile radius of UH Manoa during 90-days, spread across 12 months
  * We want to even out sampling to avoid any bias
   * Ex. We may underestimate # citations during summer or overestimate # citations during the holiday season

* What does it mean to estimate the distribution parameter based on maximum likelihood?


```python
citations_data = pd.read_csv("data/citations_counts.tsv", index_col="Day")
citations_data.head()
```




<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>Counts</th>
    </tr>
    <tr>
      <th>Day</th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>27</td>
    </tr>
    <tr>
      <th>1</th>
      <td>18</td>
    </tr>
    <tr>
      <th>2</th>
      <td>30</td>
    </tr>
    <tr>
      <th>3</th>
      <td>13</td>
    </tr>
    <tr>
      <th>4</th>
      <td>20</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize=(12,8))



plt.hist(citations_data["Counts"], bins=8, density=True, edgecolor='black', linewidth=1.2, alpha=0.5)
x = np.arange(30)



```

### Fitting the Data

* We see that the most conmmon value is $\approx$ 17
  * Poisson distribtion  with a mean ($\lambda$) of 8  is not likely to have generate this dataset
  * Poisson distribtion  with a mean ($\lambda$) of 18  is more likely to have generate this dataset that a value of $\lambda =8$
  
  
  

  





```python
plt.figure(figsize=(12,8))

from scipy.stats import poisson


np.random.seed(4)


plt.hist(citations_data["Counts"], bins=8, normed=True, edgecolor='black', linewidth=1.2)
x = np.arange(30)


p_y = poisson.pmf(x, 9)
plt.plot(x, p_y, color= 'r', linewidth=4, label="$\lambda=9$")
p_y = poisson.pmf(x, 18)
plt.plot(x, p_y, color= 'k', linewidth=4, label="$\lambda=18$")
plt.legend()
```


```python
plt.figure(figsize=(12,8))

plt.hist(citations_data["Counts"], bins=8, normed=True, edgecolor='black', linewidth=1.2)

for i in [12,13,14,15,16, 17,18,19,20]:
    p_y = poisson.pmf(x, i)
    plt.plot(x, p_y, color= 'r', linewidth=0.5, label="$\lambda=%sa$" % i, alpha=.9)

p_y = poisson.pmf(x, 17)

plt.plot(x, p_y, color= 'r', linewidth=1.5, label="$\lambda=%sa$" % i)
    
```

### Computing the Probability of Observing these Data 

* We can compute the probability given each possible parameter  
  * We limit ourselves to reasonable values
  

  * We assume the data is generated from a Poisson distribution
  * We know how to compute the probability (pmf) of a single event
          `poisson.pmf(x, lambda)`
          
  * Since all events are independent, the probability of the complete dataset is simply the product of all probabilities
              `np.product(poisson.pmf(x, lambda))`
  



### Compound probability of Independent Events

* Given a fair coin $p(H)=p(T) = 1/2$
* If we flip the coin twice, the sample sapce is \{HH, HT, TH, TT\} and  probabilities of the different outcomes is:

$$
p(HH) = p(H) * p(H) = 1/4 \\
p(HT) = p(H) * p(T) = 1/4 \\
p(TH) = p(T) * p(H) = 1/4 \\
p(TT) = p(T) * p(T) = 1/4 \\
$$

* As expected, all the probabilities sum to one 

* This only works if the probabilities are independent; i.e., obtaining `head` in the first toss does not affect the outcome of the second coin flip

### Gambler's Fallacy - Cont'd





### Gambler's Fallacy - Cont'd




### Computing the Likelihood (Probability)

* If the events are independent, the probability of the data is is the product of the probabilities

* This probability is the likelihood function $L$ of $\lambda$

$$ 
\begin{align}
L([x_1,x_2,...x_{90}] | \lambda) &= pmf(x_1, \lambda) \times  pmf(x_2, \lambda)  ....  pmf(x_{90}, \lambda) \\
  &=  \prod_{i=1}^{90} pmf(x_i, \lambda) 
\end{align}
$$




```python
data_point = 12
_lambda = 16
poisson.pmf(data_point, _lambda)
```




    0.06612887239464967




```python
data_point = 16
poisson.pmf(data_point, _lambda)
```




    0.09921753162215541




```python

poisson.pmf(citations_data["Counts"], _lambda)

```




    array([0.00335386, 0.08300552, 0.00056393, 0.08138938, 0.05591951,
           0.04959665, 0.09921753, 0.06612887, 0.09338121, 0.06989938,
           0.09921753, 0.04959665, 0.09921753, 0.09921753, 0.09921753,
           0.06989938, 0.09301644, 0.08300552, 0.09338121, 0.0309857 ,
           0.08300552, 0.09921753, 0.02155527, 0.05591951, 0.0309857 ,
           0.09301644, 0.09338121, 0.09338121, 0.09921753, 0.01198747,
           0.04959665, 0.09301644, 0.00599374, 0.08300552, 0.09338121,
           0.09921753, 0.09921753, 0.08300552, 0.09301644, 0.09301644,
           0.09921753, 0.09338121, 0.08138938, 0.09921753, 0.08300552,
           0.09301644, 0.09338121, 0.06989938, 0.09921753, 0.06612887,
           0.09921753, 0.00105737, 0.09921753, 0.09921753, 0.08300552,
           0.09921753, 0.09301644, 0.06989938, 0.06612887, 0.02131106,
           0.05591951, 0.09301644, 0.09301644, 0.09301644, 0.02131106,
           0.09338121, 0.0309857 , 0.00919692, 0.06612887, 0.08138938,
           0.0340977 , 0.01437018, 0.09338121, 0.05591951, 0.0309857 ,
           0.09921753, 0.09301644, 0.04260534, 0.06612887, 0.04260534,
           0.09921753, 0.0340977 , 0.08138938, 0.06612887, 0.00919692,
           0.00919692, 0.09338121, 0.09921753, 0.09921753, 0.06612887])




```python
np.product(poisson.pmf(citations_data["Counts"], _lambda))

```




    1.1152800100792753e-113



### Understanding the Likelihood (Probability)

1. What does the `np.product(poisson.pmf(poi_data, 16))` compute?


2. Why is the resulting probability so small?
  * Thik, for instance about a sample form 3 days (16,15,16)
  * The sample space is large and includes (0,1,8), (1,10,10), (21,11,28), ...
  
  
3. The probability of observing (12,11,11) is just one event (possibility)  among  many (countably infinite) other events




```python
_lambda = 16
np.product(poisson.pmf((16,15,16), _lambda))
```




    0.0009767091466218046




```python
np.product(poisson.pmf((9,11,11), 16))
```




    5.24215504845479e-05




```python
total_prob = 0
for i in range(0, 35):
    for j in range(0, 35):
        for k in range(0, 35):
            total_prob += np.product(poisson.pmf((i,j,k), 16))
            
print(total_prob)
```

    0.9999192928464893



```python
for i in range(12, 21):
    likelihood_data_i = np.product(poisson.pmf(citations_data["Counts"], i))
    print("The likelihood of the data given lambda=%s is %s" % (i, likelihood_data_i))
```

    The likelihood of the data given lambda=12 is 7.225758340021846e-142
    The likelihood of the data given lambda=13 is 1.3062039910366992e-129
    The likelihood of the data given lambda=14 is 3.6837061119261494e-121
    The likelihood of the data given lambda=15 is 5.438292115885742e-116
    The likelihood of the data given lambda=16 is 1.1152800100792753e-113
    The likelihood of the data given lambda=17 is 7.058672617915774e-114
    The likelihood of the data given lambda=18 is 2.670870493915664e-116
    The likelihood of the data given lambda=19 is 1.051303302700141e-120
    The likelihood of the data given lambda=20 is 6.877960001068855e-127



```python
plt.figure(figsize=(12,4))

x = range(12, 21)
l_x = [np.product(poisson.pmf(citations_data["Counts"], i)) for i in x]


frame1 = plt.gca()
frame1.axes.yaxis.set_ticklabels([])
plt.scatter(x, l_x)

_= plt.ylim(-1.1152800100792753e-113, 1.1152800100792753e-112)

```

### Understanding the Likelihood - Cont'd


* Working with small numbers can lead to likelihood to be very small
  * Can lead to arithmetic floating point underflow 
    * Results of a calculation smaller than what can actually be represented in memory
* Rather than computing likelihood, we can calculate Log-likelihood
  * Scale changes but share of the function does not
  * Those are no longer probabilities as they don't sum to 1

We know that:

$$
log(x \dot y) = log(x) + log(y)
$$

Therefore:


\begin{align}
log(L([x_1,x_2,...x_90] | \lambda)) &= log(pmf(x_1, \lambda)) \times  log(pmf(x_2, \lambda))  ....  log(pmf(x_{90}, \lambda)) \\
  &=  \sum_{i=1}^{90} log(pmf(x_i, \lambda)) 
\end{align}



```python
plt.figure(figsize=(12,8))

x = range(12, 21)
l_x = [np.sum(np.log(poisson.pmf(citations_data["Counts"], i))) for i in x]

plt.scatter(x, l_x)

```

### Maximum likelihood For the Poisson Distribution

* We hypothesize that the data was generated using a parameter $\lambda$ 
  * Among all possible $lambda$ there is at least one under which the data is most likely
  * We say that this  value maximizes the likelihood

* We can quickly show that the value that maximized the likelihood merely is the sample mean!

```python 
>>> np.mean(citations_data["Counts"]) 
16.41
```
* While we assumed that the data is best represented using a Poisson, we could have manually tested other distributions and selected the one that provides the highest likelihood
 
  * We will discuss "goodness of fit," which tackles the issue of selecting the most appropriate model for the data, later


# Example: A/B Testing


* In an A/B test, we are interested in comparing two versions of a webpage
  * We would like to determine which one performs better at a task `X`
    
* Ex.: Is Version A (new) or Version B (old) better at signing up new users?





### Maximum Likelihood For a Binomial Distribution

* For now, we are only concerned with version A.

* We test 8 individuals (unrealistically small sample) using version A
 * 5 sign up for the newsletter

*  What can we infer about the type of data and parameter that gave rise to this outcome
  
1. The data is distributed according to a binomial distribution 
   * $n=8$ and $p \in [0,1]$

2 Intuitively, we know that parameter cannot too small or too big
  * From the data, p is somewhere in the upper  side of the range
  * Ex the data is more likely with a $p=0.8$ than with  $p=0.03$

* We can determine $p$ using the same approach used with traffic citations data


```python
from scipy.stats import binom
n = 8 
p = 0.5
binom.pmf(3,n ,p)
```




    0.21875000000000014




```python
binom.pmf(4, n ,p)
```




    0.27343750000000017




```python
x = np.arange(0.001,1,0.025)
l_x_vals = []
for p in x:
    l_x_i = np.sum(np.log(binom.pmf(5, n , p)))
    l_x_vals.append(l_x_i)
    
plt.scatter(x, l_x_vals)
plt.scatter(x[np.argmax(l_x_vals)], np.max(l_x_vals))

plt.axvline(0.625, color='r')

```

### Maximum Likelihood For a Binomial Distribution - Cont'd

* The maximum likelihood is consistent with our intuition

  * If we toss a coin 100 times and observe 85 heads, then we assume the probability of heads is 0.85 
  * Since 5 out of clients signed up, we assume that $p = 5/ 8 =  0.625$ 

* This can also be easily proven analytically



###  Maximum Likelihood Estimation

* The maximum likelihood estimators for paramters of common distributions can be computed analytically

Ex. The maximum estimator of the mean in a gaussian is simply the sample mean and the estimator of the std. deviation is the sample's standard deviation
    
    
$$
\hat{\mu} = \frac{1}{n} \sum_{i=1}^{n}x_i, 
$$
and
$$
\hat{\mu} = \frac{1}{n} \sum_{i=1}^{n}{(x_i - \bar{x})^2}.
$$
    
    
* Recall that ML estimators are simply point-estimates that maximize the likelihood of data

* Do not give us a probability of the value we generate
 * We do not have the probability or confidence that it is the correct parameter of the population


### Parameter Estimation and Variance


* If we have the relevant data for the entire population, we can simply calculate the parameter. 
* We know that a statistic  based on a random sample can be a reasonable estimate of an unknown parameter in the population
 * BUT: We cannot guarantee that two samples drawn at random will have exctly the same paramter
 
 
* The questions i then: ** "How different could this estimate have been, if the sample had come out differently?"**
 

###  Sampling and Resampling

* A random sample is just one of numerous possible random samples

  * Our parameter estimate is, therefore, just one of numerous plausible estimates


* By how much could those estimates vary? 

* Good was to test this is by drawing another sample from the population and compare the new estimate to the old one

  * Sampling from the population can be costly, tedious or simply not possible



### The Bootstrap: Resampling from the Sample ###

* Key Idea: A large random sample is likely to resemble the population from which it is drawn. 

 * This observation allows data scientists replicate the sampling procedure by *sampling from the sample*. 

* Here are the steps of *the bootstrap method* for generating another random sample that resembles the population:

- **Treat the original sample as if it were the population.**

- **Draw from the sample**, at random **with** replacement, **the same number of times as the original sample size**

* For instance, rather than sample new citation, we consider our dataset as a population  and draw new samples from it


```python
test_data = [1,2,3,4]

np.random.choice(test_data, 4, replace=True)

# np.random.choice uses replace = True by default
```




    array([1, 4, 2, 3])




```python
plt.figure(figsize=(12,8))



plt.hist(citations_data["Counts"], bins=8, normed=True, edgecolor='black', linewidth=1.2, alpha =0.5)
x = np.arange(30)
bootstrap_1 = np.random.choice(citations_data["Counts"], len(citations_data["Counts"]))
_ = plt.hist(bootstrap_1, bins=8, normed=True, edgecolor='black', linewidth=1.2, alpha =0.5)


```


```python
print(np.mean(bootstrap_1))
```

    16.966666666666665



```python
[np.mean(np.random.choice(citations_data["Counts"], len(citations_data["Counts"]))) for _ in range(10)]
```




    [16.733333333333334,
     17.011111111111113,
     16.511111111111113,
     15.911111111111111,
     16.11111111111111,
     16.322222222222223,
     16.58888888888889,
     16.377777777777776,
     16.844444444444445,
     17.933333333333334]



### Remember 
1 - Resample the same number of times as the original sample size.  

  - The reason is that the variability of an estimate depends on the size of the sample. 
  - To see how different the sample could have been, we have to compare it to the median of other samples of size 500.

- If we drew 500 times at random *without* replacement from our sample of size 500, we would just get the same sample back. 

  - By drawing *with* replacement, we create the possibility for the new samples to be different from the original, because some employees might be drawn more than once and others not at all.


### Why is This a Good Idea? 

* Similarity by transitivity

  *  By the law of averages, the distribution of the original sample is likely to resemble the population, 

  *  The distributions of all the "resamples" are likely to resemble the original sample.

  * So the distributions of all the resamples are likely to resemble the population as well. 

### Resampling 






```python

```

15c941eb71ac837de2f2289f056f5abad1c6cd13
ipynb
Jupyter Notebook
morea/parameter_estimation/resources/.ipynb_checkpoints/param_esitmation_maximum_likelihood-checkpoint.ipynb
mahdi-b/ICS434-Sp2022
b2a6daf9833f4f4ec67b74c81794b1d7259b9909
Apache-2.0
morea/parameter_estimation/resources/.ipynb_checkpoints/param_esitmation_maximum_likelihood-checkpoint.ipynb
mahdi-b/ICS434-Sp2022
b2a6daf9833f4f4ec67b74c81794b1d7259b9909
Apache-2.0
morea/parameter_estimation/resources/.ipynb_checkpoints/param_esitmation_maximum_likelihood-checkpoint.ipynb
mahdi-b/ICS434-Sp2022
b2a6daf9833f4f4ec67b74c81794b1d7259b9909
Apache-2.0
2022-03-26T04:28:08.000Z
2022-03-26T04:28:08.000Z
Qwen/Qwen-72B
1. YES
2. YES