# DATA 5600: Introduction to Regression and Machine Learning for Analytics

## __Computational Thinking, Statistical Thinking, and First Steps with the Monte Carlo Method__ <br>

Author:      Tyler J. Brough <br>
Updated: September 22, 2021 <br>

---

<br>


```python
import numpy as np
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt
from typing import List

sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = [10, 5]
```


```python
students = ['Tyler Clayson', 'Marshall Day', 'Liz Giles', 'Connor Hoiland', 'Josie Hull',
            'Brian Kissmer', 'Andrew Martin', 'Jacob Needham', 'Tori Rallison', 'Emily Rice',
            'Chris Simiskey', 'Wes Smith', 'Andrew Walker', 'Austin Warr', 'Connor Waterman']


def draw_a_student(students: List[str]): 
    the_student = np.random.choice(students)
    print(f"{the_student}, come on down!")
```


```python
draw_a_student(students)
```

    Jacob Needham, come on down!


<br>

---

### __Example 1:__ The Binomial Distribution

#### __Preliminary Bits__

<br>

__Bernoulli Trials:__ A set of $n$ independent binary variables in which the $j$th observation is either a 'success' or 'failure', with probability of success, $\theta$, being the same for all trials.

<br>

#### The __Binomial distribution__


<br>

The distribution of the number of _successes_, $X$, in a series of $n$-independent ___Bernoulli trials___ where the probability of success at each trial is $\theta$ and the probability of failure is $1 - \theta$. Specifically the distribution is given by

<br>

$$
\large{Pr(X = x) = \frac{n!}{x!(n-x)!} \theta^{x}(1 - \theta)^{n-x}, \quad x = 1, 2, \ldots, n}
$$

<br>

The mean, variance, skewness and kurtosis of the distribution are as follows:

<br>

$$
\begin{align}
\mbox{mean} &= n\theta \\
& \\
\mbox{variance} &= n\theta(1 - \theta) \\
& \\
\mbox{skewness} &= ((1 - \theta) - \theta) / (n\theta(1-\theta))^{\frac{1}{2}} \\
& \\
\mbox{kurtosis} &= 3 - \frac{6}{n} + \frac{1}{n\theta(1 - \theta)} \\
\end{align}
$$

<br>


```python
np.random.binomial(n=10, p=0.5, size=1)
```




    array([3])




```python
M = 10_000
N = 50
theta = 0.5

heads = np.zeros(M)

## the main simulation loop
for i in range(M):
    heads[i] = np.random.binomial(n=N, p=theta)
    
(np.mean(heads), np.std(heads, ddof=1))
```




    (24.9841, 3.585238148750837)




```python
N * theta
```




    25.0




```python
np.sqrt(N * theta * (1. - theta))
```




    3.5355339059327378




```python
plt.hist(heads, bins=100);
```

<br>

__NB:__ We can actually vectorize this simulation quite simply as the following:

<br>


```python
more_heads = np.random.binomial(n=N, p=theta, size=M)
```


```python
plt.hist(more_heads, bins=25);
```

<br>

__Q:__ What's the difference?

<br>


##### __Standarizing the Results__

<br>


```python
z_scores = (more_heads - np.mean(more_heads)) / np.std(more_heads, ddof=1)
plt.hist(z_scores, bins=100);
```


```python
(np.mean(z_scores), np.std(z_scores, ddof=1))
```




    (4.092726157978177e-16, 1.0)



<br>

---

<br>

### __Example 2:__ The Poisson Distribution

<br>

#### The __Poisson Distribution__ 

<br>

The probability distribution of the number of occurrences, $X$, or some random event, in an interval of time or space. Given by

<br>

$$
\large{Pr(X = x) = \frac{e^{-\lambda} \lambda^{x}}{x!}, \quad x = 0, 1, \ldots}
$$

<br>

The mean and variance of the distribution are both $\lambda$. The skewness of the distribution is $1 / \sqrt{\lambda}$ and its kurtosis is $3 + 1 / \lambda$. The distribution is positively skewed for all values of $\lambda$. 

<br>


```python
np.random.poisson?
```


    [0;31mDocstring:[0m
    poisson(lam=1.0, size=None)
    
    Draw samples from a Poisson distribution.
    
    The Poisson distribution is the limit of the binomial distribution
    for large N.
    
    .. note::
        New code should use the ``poisson`` method of a ``default_rng()``
        instance instead; please see the :ref:`random-quick-start`.
    
    Parameters
    ----------
    lam : float or array_like of floats
        Expectation of interval, must be >= 0. A sequence of expectation
        intervals must be broadcastable over the requested size.
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  If size is ``None`` (default),
        a single value is returned if ``lam`` is a scalar. Otherwise,
        ``np.array(lam).size`` samples are drawn.
    
    Returns
    -------
    out : ndarray or scalar
        Drawn samples from the parameterized Poisson distribution.
    
    See Also
    --------
    Generator.poisson: which should be used for new code.
    
    Notes
    -----
    The Poisson distribution
    
    .. math:: f(k; \lambda)=\frac{\lambda^k e^{-\lambda}}{k!}
    
    For events with an expected separation :math:`\lambda` the Poisson
    distribution :math:`f(k; \lambda)` describes the probability of
    :math:`k` events occurring within the observed
    interval :math:`\lambda`.
    
    Because the output is limited to the range of the C int64 type, a
    ValueError is raised when `lam` is within 10 sigma of the maximum
    representable value.
    
    References
    ----------
    .. [1] Weisstein, Eric W. "Poisson Distribution."
           From MathWorld--A Wolfram Web Resource.
           http://mathworld.wolfram.com/PoissonDistribution.html
    .. [2] Wikipedia, "Poisson distribution",
           https://en.wikipedia.org/wiki/Poisson_distribution
    
    Examples
    --------
    Draw samples from the distribution:
    
    >>> import numpy as np
    >>> s = np.random.poisson(5, 10000)
    
    Display histogram of the sample:
    
    >>> import matplotlib.pyplot as plt
    >>> count, bins, ignored = plt.hist(s, 14, density=True)
    >>> plt.show()
    
    Draw each 100 values for lambda 100 and 500:
    
    >>> s = np.random.poisson(lam=(100., 500.), size=(100, 2))
    [0;31mType:[0m      builtin_function_or_method




```python
np.random.poisson(lam=4.5, size=10)
```




    array([7, 6, 7, 5, 5, 4, 5, 9, 1, 4])




```python
draw_a_student(students)
```

    Liz Giles, come on down!



```python
M = 10_000
N = 100
lam = 2

avg_counts = np.zeros(M)
for i in range(M):
    counts = np.random.poisson(lam=lam, size = N)
    avg_counts[i] = np.mean(counts)
    
```


```python
avg_counts[:10]
```




    array([1.81, 1.88, 1.87, 1.98, 2.02, 2.02, 1.93, 2.08, 1.81, 2.05])



<br>

---

<br>


```python
plt.hist(avg_counts, bins=25);
```


```python
z_scores = (avg_counts-lam)/np.sqrt(lam)
plt.hist(z_scores, bins=50);
```

### __Example 3:__ The Exponential Distribution

<br>


```python

```

<br>

#### __Preliminary Bits__

<br>

__Poisson Process:__ A point process with independent increments at constant intensity, say $\lambda$. The count after time $t$ has therefore a _Poisson distribution_ with mean $\lambda t$. 

<br>

__Point Process:__ A stochastic process concerned with the occurrence of events at points of time or space determined by some chance mechanism. 

<br>

#### The __Exponential Distribution__

<br>

The probability distribution $f(x)$ given by

<br>

$$
\large{f(x) = \lambda e^{-\lambda x}, \quad x > 0}
$$

<br>


The mean, variance, skewness, and kurtosis of the distribution are as follows:

<br>

$$
\begin{align}
\mbox{mean} &= 1 / \lambda \\
& \\
\mbox{variance} &= 1 / \lambda^{2} \\
& \\
\mbox{skewness} &= 2 \\
& \\
\mbox{kurtosis} &= 9 \\
& \\
\end{align}
$$

<br>

__NB:__ The distribution of intervals between independent consecutive random events, i.e. those following a ___Poisson process___.

<br>

<br>

__NB:__ `numpy` uses a different parameterization from the one above:

<br>

$$
f(x) = \frac{1}{\beta} e^{-\frac{1}{\beta}x}, \quad \mbox{which just substitutes} \quad \lambda = \frac{1}{\beta}
$$

<br>


```python
np.random.exponential?
```


```python
np.random.exponential(scale = 4, size=10)
```


```python
draw_a_student(students)
```


```python

```


```python

```


```python

```

<br>

---

<br>

### __A Simple Linear Regression: The Slope Coefficient__

<br>

Consider the following simple linear regression:

<br>

$$
\large{y_{i} = a + b x_{i} + \varepsilon_{t}}
$$

<br>

We are interested in the OLS estimate of the $b$ coefficient, $\hat{b}$, and its sampling distribution.

<br>


```python

```

2700e185134c0a16f0320e736423248f079a066c
ipynb
Jupyter Notebook
Lectures/Gelman et al/Chapter04/Computational-Thinking.ipynb
CAIJedi/DATA-5600-FL21
df88e9462eb87c277e52966359638f07874b60fc
CC0-1.0
2021-11-11T00:21:39.000Z
2021-11-11T00:21:39.000Z
Lectures/Gelman et al/Chapter04/Computational-Thinking.ipynb
CAIJedi/DATA-5600-FL21
df88e9462eb87c277e52966359638f07874b60fc
CC0-1.0
Lectures/Gelman et al/Chapter04/Computational-Thinking.ipynb
CAIJedi/DATA-5600-FL21
df88e9462eb87c277e52966359638f07874b60fc
CC0-1.0
2021-09-15T19:40:06.000Z
2021-11-09T02:48:36.000Z
Qwen/Qwen-72B
1. YES
2. YES