import numpy as np

# Data
age = np.array([13, 14, 14,12, 9, 15, 10, 14, 9, 14, 13, 12, 9, 10, 15, 11, 15, 
                11, 7, 13, 13, 10, 9, 6, 11, 15, 13, 10, 9, 9, 15, 14, 14, 10, 14, 11, 13, 14, 10])
price = np.array([2950, 2300, 3900, 2800, 5000, 2999, 3950, 2995, 4500, 2800, 1990, 
                  3500, 5100, 3900, 2900, 4950, 2000, 3400, 8999, 4000, 2950, 3250, 
                  3950, 4600, 4500, 1600, 3900, 4200, 6500, 3500, 2999, 2600, 3250, 
                  2500, 2400, 3990, 4600, 450,4700])/1000.


from pymc import Normal, Gamma, deterministic, MCMC, Matplot

# Constant priors for parameters
a = Normal('a', 0, 0.0001)
b = Normal('b', 0, 0.0001)

# Precision of normal distribution of prices
tau = Gamma('tau', alpha=0.1, beta=0.1)

@deterministic
def mu(x=age, a=a, b=b):
    # Linear age-price model
    return a + b*x

# Sampling distribution of prices
p = Normal('p', mu, tau, value=price, observed=True)

M = MCMC(locals(), db='pickle')
M.sample(iter=20000, burn=10000)

%matplotlib inline
Matplot.plot(b)

import matplotlib.pyplot as plt
from pymc.examples.disaster_model import disasters_array

plt.figure(figsize=(12.5, 3.5))
n_count_data = len(disasters_array)
plt.bar(np.arange(1851, 1962), disasters_array, color="#348ABD")
plt.xlabel("Year")
plt.ylabel("Disasters")
plt.title("UK coal mining disasters, 1851-1962")
plt.xlim(1851, 1962);

from pymc import DiscreteUniform, Exponential, Poisson, deterministic

switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110)

switchpoint

early_mean = Exponential('early_mean', beta=1., value=3)
late_mean = Exponential('late_mean', beta=1., value=1)

@deterministic
def rate(s=switchpoint, e=early_mean, l=late_mean):
    # Create a vector of Poisson means
    out = np.empty(len(disasters_array))
    out[:s] = e
    out[s:] = l
    return out

disasters = Poisson('disasters', mu=rate, value=disasters_array, 
                    observed=True)

switchpoint.parents

disasters.parents

rate.parents

rate.children

from pymc import graph, MCMC
graph.dag(MCMC([switchpoint, rate, early_mean, late_mean, disasters]))

!dot MCMC.dot -Tpng -o images/dag.png
from IPython.core.display import Image
Image('images/dag.png')

disasters.value

switchpoint.value

early_mean.value

late_mean.value

rate.value

switchpoint.logp

late_mean.value = 5

disasters.logp

early_mean.logp

late_mean.logp

from pymc.examples import disaster_model
M = MCMC(disaster_model)

M.sample(iter=10000, burn=1000)

M.trace('switchpoint')[1000:]

M.early_mean.trace()

plt.hist(M.trace('late_mean')[:])

from pymc.Matplot import plot
plot(M)

Matplot.summary_plot([M.early_mean, M.late_mean])

M.early_mean.stats()

x = np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, None, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, None, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])

masked_values = np.ma.masked_values(x, value=None)
masked_values

disasters = Poisson('disasters', mu=rate, 
                    value=masked_values, observed=True)

def missing_data_model():

    # Switchpoint
    switch = DiscreteUniform('switch', lower=0, upper=110)
    # Early mean
    early_mean = Exponential('early_mean', beta=1)
    # Late mean
    late_mean = Exponential('late_mean', beta=1)
    
    @deterministic(plot=False)
    def rate(s=switch, e=early_mean, l=late_mean):
        """Allocate appropriate mean to time series"""
        out = np.empty(len(disasters_array))
        # Early mean prior to switchpoint
        out[:s] = e
        # Late mean following switchpoint
        out[s:] = l
        return out
    
    masked_values = np.ma.masked_values(x, value=None)
    
    # Pass masked array to data stochastic, and it does the right thing
    disasters = Poisson('disasters', mu=rate, value=masked_values, observed=True)
    
    return locals()

M.early_mean.summary()

M_missing = MCMC(missing_data_model())
M_missing.sample(5000)

Matplot.plot(M_missing.disasters)

M.step_method_dict

from pymc import Slicer
M.use_step_method(Slicer, disaster_model.late_mean, w=0.5)

from IPython.core.display import HTML
def css_styling():
    styles = open("styles/custom.css", "r").read()
    return HTML(styles)
css_styling()