<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The Dirichlet process mixture model is incredibly flexible in terms of the family of parametric component distributions {fθ | fθ∈Θ}{fθ | fθ∈Θ}. We illustrate this flexibility below by using Poisson component distributions to estimate the density of sunspots per year.
Step1:  Generating data
Step2:  Model specification
Our initial beliefs about the parameters are quite informative (sd=1) and a bit off the true values.
We will use the model
Step3:  For the sunspot model, the posterior distribution of αα is concentrated between 0.6 and 1.2, indicating that we should expect more components to contribute non-negligible amounts to the mixture than for the Old Faithful waiting time model.
Step4:  Indeed, we see that between ten and fifteen mixture components have appreciable posterior expected weight.
Step5:  We now calculate and plot the fitted density estimate.
Step6:  Again, we can decompose the posterior expected density into weighted mixture densities.
<ASSISTANT_TASK:>
Python Code:
# pymc3.distributions.DensityDist?
import matplotlib.pyplot as plt
import matplotlib as mpl
from pymc3 import Model, Normal, Slice
from pymc3 import sample
from pymc3 import traceplot
from pymc3.distributions import Interpolated
from theano import as_op
import theano.tensor as tt
import numpy as np
from scipy import stats
from matplotlib import pyplot as plt
import numpy as np
import pymc3 as pm
import scipy as sp
import seaborn as sns
from statsmodels.datasets import get_rdataset
from theano import tensor as tt
%matplotlib inline
%load_ext version_information
%version_information pymc3, statsmodels, pandas
Explanation: The Dirichlet process mixture model is incredibly flexible in terms of the family of parametric component distributions {fθ | fθ∈Θ}{fθ | fθ∈Θ}. We illustrate this flexibility below by using Poisson component distributions to estimate the density of sunspots per year.
End of explanation
sunspot_df = get_rdataset('sunspot.year', cache=True).data
sunspot_df.head()
sunspot_df.plot(x='time')
Explanation: Generating data
End of explanation
SEED = 8675309 # from random.org
np.random.seed(SEED)
K = 50
N = sunspot_df.shape[0]
def stick_breaking(beta):
    portion_remaining = tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]])
    return beta * portion_remaining
with pm.Model() as model:
    alpha = pm.Gamma('alpha', 1., 1.)
    beta = pm.Beta('beta', 1, alpha, shape=K)
    w = pm.Deterministic('w', stick_breaking(beta))
    mu = pm.Uniform('mu', 0., 300., shape=K)
    obs = pm.Mixture('obs', w, pm.Poisson.dist(mu), observed=sunspot_df['sunspot.year'])
with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, step=step, tune=90000, random_seed=SEED, njobs=6)
Explanation: Model specification
Our initial beliefs about the parameters are quite informative (sd=1) and a bit off the true values.
We will use the model:
$\alpha \sim Gamma(1,1) \
\beta_1,...,\beta_k \sim Beat(1,\alpha) \
w_i = \beta_i \prod^i_{j=i-1}(1-\beta_j \
\lambda_i,...,\lambda_k \sim U(0,300) \
x | w_i,\lambda_i \sim \sum^K_{i=1}w_i Poisson(\lambda_i)$
End of explanation
pm.traceplot(trace, varnames=['alpha']);
Explanation: For the sunspot model, the posterior distribution of αα is concentrated between 0.6 and 1.2, indicating that we should expect more components to contribute non-negligible amounts to the mixture than for the Old Faithful waiting time model.
End of explanation
fig, ax = plt.subplots(figsize=(8, 6))
plot_w = np.arange(K) + 1
ax.bar(plot_w - 0.5, trace['w'].mean(axis=0), width=1., lw=0);
ax.set_xlim(0.5, K);
ax.set_xlabel('Component');
ax.set_ylabel('Posterior expected mixture weight');
Explanation: Indeed, we see that between ten and fifteen mixture components have appreciable posterior expected weight.
End of explanation
x_plot = np.arange(250)
post_pmf_contribs = sp.stats.poisson.pmf(np.atleast_3d(x_plot),
                                         trace['mu'][:, np.newaxis, :])
post_pmfs = (trace['w'][:, np.newaxis, :] * post_pmf_contribs).sum(axis=-1)
post_pmf_low, post_pmf_high = np.percentile(post_pmfs, [2.5, 97.5], axis=0)
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist(sunspot_df['sunspot.year'].values, bins=40, normed=True, lw=0, alpha=0.75);
ax.fill_between(x_plot, post_pmf_low, post_pmf_high,
                 color='gray', alpha=0.45)
ax.plot(x_plot, post_pmfs[0],
        c='gray', label='Posterior sample densities');
ax.plot(x_plot, post_pmfs[::200].T, c='gray');
ax.plot(x_plot, post_pmfs.mean(axis=0),
        c='k', label='Posterior expected density');
ax.set_xlabel('Yearly sunspot count');
ax.set_yticklabels([]);
ax.legend(loc=1);
Explanation: We now calculate and plot the fitted density estimate.
End of explanation
fig, ax = plt.subplots(figsize=(8, 6))
ax.hist(sunspot_df['sunspot.year'].values, bins=40, normed=True, lw=0, alpha=0.75);
ax.plot(x_plot, post_pmfs.mean(axis=0),
        c='k', label='Posterior expected density');
ax.plot(x_plot, (trace['w'][:, np.newaxis, :] * post_pmf_contribs).mean(axis=0)[:, 0],
        '--', c='k', label='Posterior expected\nmixture components\n(weighted)');
ax.plot(x_plot, (trace['w'][:, np.newaxis, :] * post_pmf_contribs).mean(axis=0),
        '--', c='k');
ax.set_xlabel('Yearly sunspot count');
ax.set_yticklabels([]);
ax.legend(loc=1);
Explanation: Again, we can decompose the posterior expected density into weighted mixture densities.
End of explanation
<END_TASK>