# Acknowledgements 

This is a summary notebook borrowing in its entirety from William Koehrsen's https://williamkoehrsen.medium.com/ article in Towards Data Science and its corresponding github repo
1. https://towardsdatascience.com/the-poisson-distribution-and-poisson-process-explained-4e2cb17d459
2. https://github.com/WillKoehrsen/Data-Analysis/blob/master/poisson/poisson.ipynb

<b>A Poisson process is </b>

1. Model for discrete events
2. Average time between events is known (and constant)
3. Exact timing of events is random
4. Events are independent
5. Two events cannot occur at same time

Poisson Distribution gives the probability of a number of events in an interval generated by a Poisson process.


We use meteor observation to explain this. Some technical terms:
1. Asteroids are large chunks of rock orbiting the sun in the asteroid belt. 
2. Pieces of asteroids that break off become meteoroids. A meteoroid can come from an asteroid, a comet, or a piece of a planet and is usually millimeters in diameter but can be up to a kilometer. 
3. If the meteoroid survives its trip through the atmosphere and impacts Earth, it’s called a meteorite.
4. Meteors are the streaks of light you see in the sky that are caused by pieces of debris called meteoroids burning up in the atmosphere


<b>We can use Poisson to describe two things:</b>
1. Use Poisson Process for probablity of events in an interval
2. Find waiting time between two events using Poisson distribution

Some processes that are close approximation of Poisson Processes:
1. Meteors
2. Failures in a website
3. Customers calling for support
4. Stock price movements (we know average stock price movements per day)

<b>Two observations from the definition of Poisson process:</b>
1. Two events cannot occur at same time implies that we can think of each sub-interval of a Poisson process as a Bernoulli Trial, that is, either a success or a failure. With our website, the entire interval may be 600 days, but each sub-interval — one day — our website either goes down or it doesn’t.
2. Average time between events but they are randomly spaced (stochastic). We might have back-to-back failures, but we could also go years between failures due to the randomness of the process

\begin{equation}
\begin{aligned}
P(X=k) = \dfrac{\lambda^{k} }{k!} \, e^{-\lambda}
\end{aligned}
\tag{Equation 1}
\end{equation}

where parameter $ \lambda $ > 0. 


Let r be number of events per unit time $ \dfrac{events}{time} $. Then r can be called rate parameter. Then we have $ \lambda \, = \, rt $. Substituting  $ \lambda \, = \, rt $ in Equation 1 we get the following intuitively:


\begin{equation}
\begin{aligned}
P(k \, events \, in \, interval \, t) = \dfrac{(\dfrac{events}{time}\, * \, interval \, t)^{k} }{k!} \
e^{-(\dfrac{events}{time}\, * \, interval \, t)}
\end{aligned}
\tag{Equation 2}
\end{equation}

where $ \dfrac{events}{time}\, * \, inteval \, t $ is just called the parameter for Poisson Process. We do not use the phrase time period but use interval because the process can be applied to many things such as a  area volume based on the application of Poisson Process. The parameter $ \lambda $ can be thought as expected number of events in a interval.

Equation 2 can be formally written now as:

\begin{equation}
\begin{aligned}
P(k \, events \, in \, interval \, t) = \dfrac{(rt)^{k} }{k!} \
e^{-rt}
\end{aligned}
\tag{Equation 3}
\end{equation}


We want to plot Probability Mass Function (PMF) against events for different rate parameters. 

We have to first prepare required Python code. Here we go



```python
import pandas as pd
import numpy as np

from scipy.special import factorial # Poisson distribution formula has factorial term

# Display all cell outputs
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'

# Visualizations
import chart_studio
chart_studio.tools.set_credentials_file(username='datavictor', api_key='NZg3CXiXgzQM1e7ytw8U')
chart_studio.tools.set_config_file(world_readable=True,sharing='public')

from chart_studio import plotly # replaces import plotly.plotly as py
import plotly.graph_objs as go
from plotly.offline import iplot

# Cufflinks for dataframes
import cufflinks as cf
cf.go_offline()
cf.set_config_file(world_readable=True, theme='pearl')
```







```python
# np.random.seed(42) legacy method. do not use
rnd = np.random.RandomState(42)
rnd # rnd will now replace all access to np.random
```




    RandomState(MT19937) at 0x7F13463C5140



# 1. Calculating Possion Probability


```python
# Variables defiend in this cell are used throughout the notebook
events_per_minute = 1/12
minutes = 60
events_per_hour = events_per_minute * minutes
lambda_param = events_per_hour
print(f'Rate parameter is {lambda_param}')
```

    Rate parameter is 5.0



```python
# Calculate probability of k events in specified number of minutes
def poisson_probability(events_per_minute, minutes, k):
    lambda_param = events_per_minute * minutes
    return np.exp(-lambda_param) * np.power(lambda_param, k) / factorial(k)
```


```python
#Sample
k = 3
k_meteor_probability = poisson_probability(events_per_minute, minutes, k)
print(f'The probability of {k} meteors in {minutes} minutes is {100*k_meteor_probability:.2f}%.')
```

    The probability of 3 meteors in 60 minutes is 14.04%.



```python
# We can pass a list or numpy array as the 3rd param and we will get back a list or numpy array of probabilities
events_12 = np.arange(12)
print(type(events_12))
prob_12 = poisson_probability(events_per_minute, minutes, events_12)
type(prob_12)
print(f'The most likely value is {np.argmax(prob_12)} with probability {np.max(prob_12):.4f}')
```

    <class 'numpy.ndarray'>





    numpy.ndarray



    The most likely value is 4 with probability 0.1755


We can use the built in poisson function


```python
x = rnd.poisson(lambda_rate_param, 10000) #rnd replaces the previous np.random
print(type(x))
(x == 3).mean()
```

    <class 'numpy.ndarray'>





    0.1377



# 2. Plotting Poisson Mass Distribution

This is PMF because events are discrete


```python
# Plot PMF of Poisson distribution
def plot_pmf(x, p_x, title=''):
    df = pd.DataFrame({'x': x, 'y': p_x})
    # print(f'The most likely value is {np.argmax(p_x)} with probability {np.max(p_x):.4f}')
    annotations = [dict(x=x, y=y+0.01, text=f'{y:.2f}', 
                        showarrow=False, textangle=0) for x, y in zip(df['x'], df['y'])]
    df.iplot(kind='scatter', mode='markers+lines',
             x='x', y='y', xTitle='Number of Events',
             yTitle='Probability', annotations=annotations,
             title=title)
```


```python
plot_pmf(events_12, prob_12, title='Probability of Number of Meteors in One Hour')
```


<div>                            <div id="002b6d86-6412-4c8e-baec-aa9e9e4bd714" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>


# 3. Plotting Probability for different rates and different time period


NOTE: The most likely number of events in the interval for each curve is the rate parameter


```python
def plot_different_rates(events_per_minute, minutes, ns, title=''):
    df = pd.DataFrame()
    annotations=[]
    colors = ['orange', 'green', 'red', 'blue', 'purple', 'brown']
    for i, events in enumerate(events_per_minute):
        probs = calc_prob(events, minutes, ns)
        annotations.append(dict(x=np.argmax(probs)+1, y=np.max(probs)+0.025, 
                                text=f'{int(events * minutes)} MPH<br>Meteors = {np.argmax(probs) + 1}<br>P = {np.max(probs):.2f}',
                                color=colors[i],
                               showarrow=False, textangle=0))
        df[f'Meteors per Hour = {int(events * minutes)}'] = probs
    df.index = ns
    df.iplot(kind='scatter', mode='markers+lines', colors=colors, size=8, annotations=annotations,
             xTitle='Events', yTitle='Probability', title=title)
    return df
```


```python
df = plot_different_rates(events_per_minute=np.array([1/5, 1/12, 1/10, 1/15, 1/20, 1/30]),
                          minutes=60,
                          ns=list(range(15)), 
                          title='Probability of Meteors in 1 Hour at Different Rates')
```


<div>                            <div id="caec8002-76b8-49b7-b4c4-de1fe75fcd67" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>



```python
def plot_different_times(events_per_minute, minutes, ns, title=''):
    df = pd.DataFrame()
    annotations = []
    colors = ['orange', 'green', 'red', 'blue', 'purple', 'brown']
    for i, minute in enumerate(minutes):
        probs =  calc_prob(events_per_minute, minute, ns)
        annotations.append(dict(x=np.argmax(probs), y=np.max(probs)+0.025, 
                                color=colors[i],
                                text=f'{minute} Minutes<br>Meteors = {np.argmax(probs)}<br>P = {np.max(probs):.2f}',
                               showarrow=False, textangle=0))
        df[f'Minutes = {minute}'] = probs
    df.index = ns
    df.iplot(kind='scatter', mode='markers+lines', colors=colors,
             size=8, annotations=annotations,
             xTitle='Events', yTitle='Probability', title=title)
    return df
```


```python
df = plot_different_times(events_per_minute=1/12, minutes=np.array([30, 60, 90, 120]),
                         ns=list(range(15)), title='Probability of Meteors in Time Intervals for fixed rate of 5 meteors per hour')
```


<div>                            <div id="ae73063d-3b7c-4d64-9dde-8c51878994f6" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>


# 4. Simulation of Observations
We can use np.random.poisson to simulate 10,000 hours of observation and then make a histogram of observations. We expect to see a peak at 4 or 5 meteors since that is the most likely value.


```python
def plot_hist(x, title='',summary=True):
    df = pd.DataFrame(x)
    df.iplot(kind='hist', xTitle='Events', 
             yTitle='Count', title=title)
    if summary:
        print(df.describe())
```


```python
N = 10000
counts = np.random.poisson(lambda_rate_param, size=N)
plot_hist(counts, title=f'Distribution of Number of Meteors in 1 Hour Simulated {N} Times')
```


<div>                            <div id="047b1c88-c381-4f9e-9e0b-299bd8434d31" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>


                     0
    count  10000.00000
    mean       5.04460
    std        2.27049
    min        0.00000
    25%        3.00000
    50%        5.00000
    75%        6.00000
    max       15.00000



```python
counts = np.random.poisson(lambda_rate_param * 3, size=N)
plot_hist(counts, title=f'Distribution of Number of Meteors in 3 Hours Simulated {N} Times')
```


<div>                            <div id="097daa44-7a7e-448f-b62f-35fe12358923" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>


                      0
    count  10000.000000
    mean      14.974000
    std        3.874923
    min        3.000000
    25%       12.000000
    50%       15.000000
    75%       17.000000
    max       30.000000


# 5. Probability of Different Numbers of Events
Now let's take a look at the probability of seeing different numbers of meteors. We can find the probability by summing up the probabilities of more than a given number of events or less than or equal to a given number of events.


```python

def pr_less_than_or_equal(events_per_minute, minutes, n_query, quiet=False):
    p_n = poisson_probability(events_per_minute, minutes, np.arange(100))
    p = p_n[:n_query+1].sum() / p_n.sum()
    if not quiet:
        print(f'{int(events_per_minute*60)} Meteors Per Hour. Probability of {n_query} or fewer meteors in {int(minutes/60)} hour: {100*p:.2f}%.')
    return p

def pr_greater_than(events_per_minute, minutes, n_query, quiet=False):
    p = 1 - pr_less_than_or_equal(events_per_minute, minutes, n_query)
    if not quiet:
        print(f'{int(events_per_minute*60)} Meteors Per Hour. Probability of more than {n_query} meteors in {int(minutes/60)} hour: {100*p:.2f}%.')
    return p

assert pr_less_than_or_equal(events_per_minute, minutes, 6, True) + pr_greater_than(events_per_minute, minutes, 6, True) == 1
assert pr_less_than_or_equal(events_per_minute, minutes, 8, True) + pr_greater_than(events_per_minute, minutes, 8, True) == 1
```

    5 Meteors Per Hour. Probability of 6 or fewer meteors in 1 hour: 76.22%.
    5 Meteors Per Hour. Probability of 8 or fewer meteors in 1 hour: 93.19%.



```python
_ = pr_greater_than(events_per_minute=1/12, minutes=60, n_query=10)
```

    5 Meteors Per Hour. Probability of 10 or fewer meteors in 1 hour: 98.63%.
    5 Meteors Per Hour. Probability of more than 10 meteors in 1 hour: 1.37%.


# 6. Waiting time between Events

We would like to calculate how long we have to wait until an event occurs (This is very different from average time between two events). Derivation as follows

Reproducing the equation 3 below for reference

\begin{equation}
\begin{aligned}
P(X=k) = \dfrac{(rt)^{k} }{k!} \, e^{-rt}
\end{aligned}
\tag{Equation 3}
\end{equation}

where k is the number of events that occured between time 0 till t at the rate $ \lambda $


Let T be the time when first meteor sighting event occurs. Hence $ P(t<T) = 0 $

Similar probability can be obtained from equation 3 by substituting k = 0


\begin{equation}
\begin{aligned}
P(k=0) = \dfrac{(rt)^{0} }{0!} \, e^{-rt} \, = \, e^{-rt}
\end{aligned}
\tag{Equation 4}
\end{equation}

As an example, if the average time between meteor sightings is 12 minute, then there is a 60% chance of seeing one with wait > 6 minutes


```python

def waiting_time_more_than(events_per_minute, t, quiet=False):
    p = np.exp(-events_per_minute * t)
    if not quiet:
        print(f'{int(events_per_minute*60)} Meteors per hour. Probability of waiting more than {t} minutes: {100*p:.2f}%.')
    return p
    
def waiting_time_less_than_or_equal(events_per_minute, t, quiet=False):
    p = 1 - waiting_time_more_than(events_per_minute, t, quiet=quiet)
    if not quiet:
        print(f'{int(events_per_minute*60)} Meteors per hour. Probability of waiting at most {t} minutes: {100*p:.2f}%.')
    return p

def waiting_time_between(events_per_minute, t1, t2):
    p1 = waiting_time_less_than_or_equal(events_per_minute, t1, True)
    p2 = waiting_time_less_than_or_equal(events_per_minute, t2, True)
    p = p2-p1
    print(f'Probability of waiting between {t1} and {t2} minutes: {100*p:.2f}%.')
    return p

assert waiting_time_more_than(events_per_minute, 15, True) + waiting_time_less_than_or_equal(events_per_minute, 15, True) == 1
```


```python
_ = waiting_time_less_than_or_equal(events_per_minute, 12)
```

    5 Meteors per hour. Probability of waiting more than 12 minutes: 36.79%.
    5 Meteors per hour. Probability of waiting at most 12 minutes: 63.21%.



```python

def plot_waiting_time(events_per_minute, ts, title=''):
    p_t = waiting_time_more_than(events_per_minute, ts, quiet=True)
    
    df = pd.DataFrame({'x': ts, 'y': p_t})
    df.iplot(kind='scatter', mode='markers+lines', size=8,
             x='x', y='y', xTitle='Waiting Time',
             yTitle='Probability', 
             title=title)
    
    return p_t
```


```python
p_t = plot_waiting_time(events_per_minute, np.arange(100), title='Probability (T > t)')
```


<div>                            <div id="c6f88501-1cc4-4a7a-8e20-31a8ee5d36f9" class="plotly-graph-div" style="height:525px; width:100%;"></div>                    </div>


f96de0457b6f4a4c4b39a77e1cc5af8519d90ab4
ipynb
Jupyter Notebook
Distributions/Posson Distribution.ipynb
datavector-io/datascience
b1de0cd1c563b3c90d3f4382f0130c77e18308f9
MIT
Distributions/Posson Distribution.ipynb
datavector-io/datascience
b1de0cd1c563b3c90d3f4382f0130c77e18308f9
MIT
Distributions/Posson Distribution.ipynb
datavector-io/datascience
b1de0cd1c563b3c90d3f4382f0130c77e18308f9
MIT
Qwen/Qwen-72B
1. YES
2. YES