# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Explore properties of the Poisson distribution

# +
import seaborn as sns
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
from scipy import stats
import spacepy.toolbox as tb

sns.set(font_scale=1.5)

# -

bins = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
with pm.Model() as model:
    p1 = pm.Poisson('p1', 10) #bins, shape=len(bins))
    p2 = pm.Poisson('p2', 15) #bins, shape=len(bins))
    p1_2 = pm.Deterministic('p1_2', p1/p2)
    p12 = pm.Deterministic('p12', p1*p2)

    trace = pm.sample(10000)

pm.traceplot(trace)

pm.summary(trace)

pm.plot_posterior(trace)

# # Combine Poissons, can we get a similar answer faster than MCMC?

bins = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
with pm.Model() as model:
    p1 = pm.Poisson('p1', 10) #bins, shape=len(bins))
    p2 = pm.Poisson('p2', 15) #bins, shape=len(bins))
    p3 = pm.Poisson('p3', 21) #bins, shape=len(bins))

    p1_2 = pm.Deterministic('p1_2', (p1+p3)/p2)
    p12 = pm.Deterministic('p12', p1*p2/p3)

    trace = pm.sample(10000)

pm.plot_posterior(trace)

p1 = stats.poisson(10)
p2 = stats.poisson(15)
p3 = stats.poisson(21)

# +
# p1_2 = pm.Deterministic('p1_2', (p1+p3)/p2)
# p12 = pm.Deterministic('p12', p1*p2/p3)
XX = np.arange(50)

plt.plot(XX, p1.pmf(XX))
plt.plot(XX, p2.pmf(XX))
plt.plot(XX, p3.pmf(XX))

# -

# # How about summing a bunch of poissons

# +
np.random.seed(8675309)

vals = np.random.poisson(10, size=(16*5*3))


bins = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
with pm.Model() as model:
    mean = pm.Uniform('mean', 0, 100)
    p_dat = pm.Poisson('p_dat', mean, shape=len(vals), observed=vals)
    
    pmean = pm.Deterministic('pmean', p_dat.mean())

    trace = pm.sample(10000)
# -

pm.plot_posterior(trace)







# ## Fit a bunch of measurements to a Poisson

# +
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.special import factorial


# get poisson deviated random numbers
data = np.random.poisson(2, 1000)

# the bins should be of integer width, because poisson is an integer distribution
entries, bin_edges, patches = plt.hist(data, bins=11, range=[-0.5, 10.5], density=True)

# calculate binmiddles
bin_middles = 0.5*(bin_edges[1:] + bin_edges[:-1])

# poisson function, parameter lamb is the fit parameter
def poisson(k, lamb):
    return (lamb**k/factorial(k)) * np.exp(-lamb)

# fit with curve_fit
parameters, cov_matrix = curve_fit(poisson, bin_middles, entries) 

# plot poisson-deviation with fitted parameter
x_plot = np.linspace(0, 20, 1000)

plt.plot(x_plot, poisson(x_plot, *parameters), 'r-', lw=2)

# -

dat = np.random.poisson(10, size=(16*5*3))
sns.distplot(dat)

entries, bin_edges = np.histogram(dat, bins=11, range=[-0.5, 30], density=True)


# +
parameters, cov_matrix = curve_fit(poisson, tb.bin_edges_to_center(bin_edges), entries) 
x_plot = np.linspace(0, 30, 1000)

sns.distplot(dat, bin_edges)
plt.plot(x_plot, poisson(x_plot, *parameters), 'r-', lw=2)
# -

cov_matrix

dat.mean(), dat.mean()-cov_matrix[0][0]*dat.mean(), dat.mean()+cov_matrix[0][0]*dat.mean()

# +


bins = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
with pm.Model() as model:
    mean = pm.Uniform('mean', 0, 100)
    p_dat = pm.Poisson('p_dat', mean, observed=dat)
    

    trace = pm.sample(10000)
# -

pm.plot_posterior(trace)













StarcoderdataJupyterScriptsDedupFiltered