## Jupyter Snippet CB2nd 07_pymc

Jupyter Snippet CB2nd 07_pymc

# 7.7. Fitting a Bayesian model by sampling from a posterior distribution with a Markov Chain Monte Carlo method

``````import numpy as np
import pandas as pd
import pymc3 as pm
import matplotlib.pyplot as plt
%matplotlib inline
``````
``````# www.ncdc.noaa.gov/ibtracs/index.php?name=wmo-data
'cookbook-2nd-data/blob/master/'
'Allstorms.ibtracs_wmo.v03r05.csv?'
'raw=true',
delim_whitespace=False)
``````
``````cnt = df[df['Basin'] == ' NA'].groupby(
'Season')['Serial_Num'].nunique()
# The years from 1851 to 2012.
years = cnt.index
y0, y1 = years[0], years[-1]
arr = cnt.values
``````
``````# Plot the annual number of storms.
fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.plot(years, arr, '-o')
ax.set_xlim(y0, y1)
ax.set_xlabel("Year")
ax.set_ylabel("Number of storms")
``````

``````# We define our model.
with pm.Model() as model:
# We define our three variables.
switchpoint = pm.DiscreteUniform(
'switchpoint', lower=y0, upper=y1)
early_rate = pm.Exponential('early_rate', 1)
late_rate = pm.Exponential('late_rate', 1)
# The rate of the Poisson process is a piecewise
# constant function.
rate = pm.math.switch(switchpoint >= years,
early_rate, late_rate)
# The annual number of storms per year follows
# a Poisson distribution.
storms = pm.Poisson('storms', rate, observed=arr)
``````
``````with model:
trace = pm.sample(10000)
``````
``````Assigned Metropolis to switchpoint
Assigned NUTS to early_rate_log__
Assigned NUTS to late_rate_log__
100%|██████████| 10500/10500 [00:05<00:00, 1757.23it/s]
``````
``````pm.traceplot(trace)
``````

``````s = trace['switchpoint'].mean()
em = trace['early_rate'].mean()
lm = trace['late_rate'].mean()
s, em, lm
``````
``````(1930.171, 7.316, 14.085)
``````
``````fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.plot(years, arr, '-o')
ax.axvline(s, color='k', ls='--')
ax.plot([y0, s], [em, em], '-', lw=3)
ax.plot([s, y1], [lm, lm], '-', lw=3)
ax.set_xlim(y0, y1)
ax.set_xlabel("Year")
ax.set_ylabel("Number of storms")
``````