In [1]:

    

import pymc
import numpy as np
import pandas as pd
import seaborn as sns


    

In [32]:

    

disasters_array = 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, 2, 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, 1, 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])

switchpoint = pymc.DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')
early_mean = pymc.Exponential('early_mean', beta=1.)
late_mean = pymc.Exponential('late_mean', beta=1.)


    

In [19]:

    

@pymc.deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
    ''' Concatenate Poisson Means '''
    out = np.empty(len(disasters_array))
    out[:s] = e
    out[s:] = l
    return out


    

In [20]:

    

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


    

In [21]:

    

switchpoint.parents


    

    
Out[21]:





{'lower': 0, 'upper': 110}

In [22]:

    

disasters.parents


    

    
Out[22]:





{'mu': <pymc.PyMCObjects.Deterministic 'rate' at 0x7fbd66579c88>}

In [23]:

    

rate.children


    

    
Out[23]:





{<pymc.distributions.new_dist_class.<locals>.new_class 'disasters' at 0x7fbd6659d198>}

In [24]:

    

disasters.value


    

    
Out[24]:





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, 2, 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, 1, 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])

In [25]:

    

switchpoint.value


    

    
Out[25]:





array(63)

In [26]:

    

early_mean.value


    

    
Out[26]:





array(0.10579268109853043)

In [27]:

    

late_mean.value


    

    
Out[27]:





array(1.1865847805942225)

In [28]:

    

rate.value


    

    
Out[28]:





array([ 0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  0.10579268,  0.10579268,
        0.10579268,  0.10579268,  0.10579268,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478,  1.18658478,  1.18658478,  1.18658478,  1.18658478,
        1.18658478])

In [33]:

    

mymodel = [
    disasters_array,
    switchpoint,
    early_mean,
    late_mean,
    rate,
    disasters]


M = pymc.MCMC(mymodel)


    

In [34]:

    

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


    

    




 [-----------------100%-----------------] 10000 of 10000 complete in 0.6 sec

In [35]:

    

M.trace('switchpoint')[:]


    

    
Out[35]:





array([ 40,  87,  95,  65,  48,  64,  90, 109,   6, 103,  74,  31,  37,
        25,  81,  24,  31,  90,  83,  23,   0,  90,   4,  12,  27,  96,
        99,   4,  76,  27,  81, 104,  15,  51,  79,  65,  86,  58,  10,
       102,  83,  35,  89,  65,  98,  37,  97,  39,  59,  32,   7,  23,
        36,  85, 100, 109,   3,  71,  47,  80,  72,  79,  38, 107,  94,
        14,  52,  86,   6,  90,  71,   7, 103, 100,  93,  71, 105,  72,
       101,   2, 101,  32,  41,  30,  63,  63,  19, 104,   3, 106,  21,
        26, 105,   2,  53,  46,  87,  84,  59,  89,   1,  95,  20,  95,
        14, 107,  49,  73,  34,  22,  98,  96,  58,  47,  68,  41,  56,
        50, 108,  56,  13,  66,  16,  20,  74,  26, 102,  46,  60,  40,
        29,  23, 108,  68, 108,  50,  51,  64,  29,  38,  46, 103,  47,
        91,  46, 107,  71,  93,  37, 104,  42,  48,  20,  10,  55,  44,
        27,   3,  92,   6,  75,  41,  49,   4,  76,  35,  86,  61, 103,
       108, 106,  21,  83,  56, 108,  74, 108,  76,  95,  31,  60,  99,
        68,  29,  84,  17,  70,  30,   2,  39,  48,  41,  82,  68,  23,
        38, 103,  68,  77,  81,  61,   0,  89,  30,  75,  10,  23,  96,
        70,  31,  79,  35,   7,  47,  78, 103,  39,  72,  97,  67,  16,
        22,  82,  30,  27,  20,  97,  55,  52,  46,  38,  83,  24,  62,
        32,  86,  64,   1,  21,   4,  20,  23,  77,   4,  99,  71,  41,
        96,  71,  53,  15,  18,  91,  28,  88,  17,  84, 103,  66,  34,
        66,   5,  53,  42,  70,  62,  11,  43,  62,   3,  16,  11,  10,
       109,  89,  25,  38,  35,   3, 102,   6,  68,  61,  82, 100, 108,
        57,  20,  32, 100,  24,  88,  30,  69,  68,  62,  21,  76,  30,
        31,  12,  17,  38,   9,  21,  24,  90,  48,  52,   4,  64, 108,
        64,   7,   1,  54, 109,  70, 101,  49,  70,  65,  81,  75,  86,
        20,  72,  74, 106,  93,  36,  54,  40,  40,  81,  23,  31,  21,
        11,   3,  41,  39, 105,   0,  98,   2,  75,  94,  97,  77,  96,
         4,  63,  20,  30,  81, 100,   7,  63,  13,  88,   0,  66,  63,
        48,  15,  75, 104,   4,  74,  91,  50,  85,  59, 105,  73,  34,
        28, 107,  50,  97,  99, 107,  85, 109,  50,  39,  15,  59,  31,
        41,   5,  86,  23,   8,  31,  32,  71,  69,  81, 108,  51,  28,
         2,  26,  17,  62,  20,  31,  53,  57,  39,  32, 102,  11,  28,
        19, 102,  73,  11,  64, 110,  47,   9,  45,  18,   9,  29,  81,
        41,  15,  33, 100,  81,  55,  69,  95,  16,  24,  68,  94,  49,
        11, 110,  15,  94, 110,  79,  39,  42,  23,  70,  65,   9,  23,
        87,  45, 109,  98,  96,  88,  26,  15,   8, 109,  83,  54,  63,
        50,  73,  11,   2,  61,  94,  65, 107,  50,  78,  26,  42,  40,
       109,  61,  80,  80,  32,  93,  94,   3,  46,  91,   8,  34,   6,
        65,  49,  12,  31,   0,   7,  23,  66,  76,  89, 102,  69,  26,
        69,  30,  46,  35,  28,  45,  87,  16,  57,  90, 106,  86,   1,
       100, 105,  83,  72,  89,  94,  42,   1,  43,   6,  42,  44,  56,
        11,  52,  36,  14,  11, 102,  97,  54,  71,  82,   3,  63,  87,
       110, 107,   0,  79,   2,  16,  93,  42, 101,  58,  61,  72,  56,
        15,  51,  16,  65, 108,  77,   9,  24, 108,  87,   2,  41,  85,
        86,  32,  29, 105,  67,  54,  20,   9,  67, 107,  11,  96,  71,
        86,  93,  13,   5, 101,  42,  45,  60,  82,  49,  39, 103,  89,
       105,  90,  31,  60,  16,  31,  15,   2,   7,  65,   7,  86,  96,
        54,   5, 110,  49,  88,  68,  73,   8,  51,  33,  34,   3,  40,
        78,  19, 110,  71,  74,  12,   0,   6,  55,  55,  14,  94,  56,
        70,  87,  67,  15,  63, 109,  79,  73,  45,  93,  81,   9,  46,
        32,  35,  13,  99,  31,   3,   4, 107,  77,  10,  13, 102,  73,
        93,  49,   5,  18,   4,  59,  43,  73,  55,  66,  71,  22,  76,
        64,  73,  28,  64,  67,  53,  31,  68,  52,  17,  99,  92,  33,
        44,  19,  97,   6,  65,  35, 109,   3,  51,  45,  66,  58,  76,
         8,  27,  87,  57,   9,  21,  83,  65,  98,  40, 107,   0,  39,
        57,  61,  19,  33,  18,   1,  26,  86, 107,  29,  33,  80,  53,
        66,  52,   8,  64,  83,   8,   6,  59,  10,  38,  46,  94,  45,
        58,  66,  16,  43,  98,  24,  17,  36, 105,  18,  65,  19,  94,
        84, 102,  53,  62, 102,  65,  43,  67,  28,  57,  58,  92,  53,
        75, 106, 110,  51,  12,  73,  24,  27,   6,  93,  56,  68,  49,
        75,  29,  76,  10,  26,  99,  37,  33,  29,  82,  82, 101,  74,
        53,  65,  94,  76,  93,  67,  97,  71, 104,  93,  24,  89,  97,
        10,  42,  91,  63,  15,  79,  31,  47,  31,  83,  24,  11, 110,
        31,  71,  47,  91,  93,  88,  36,  36,  37, 105,  81,  83,   3,
        70,  83,  72,  42, 102,  68, 100,  97, 105, 105,  14, 108,  53,
        71,  68,  45,  28,  65,  67,  84,  88, 110,  51,  80,  84,  56,
        80, 108,   6,   5,  91,  79,  56,  90, 108,  54,   8,  71,  66,
       110,  23,  24,  60,  44,  41,  34,  53,  22,  74,  84,  88,   0,
        69,  73,  93,  37,  41,  22,  68,  90,  52,  69,  86,  34,  38,
        96,  29,  30])

In [38]:

    

import matplotlib.pyplot as plt
%matplotlib inline

sns.distplot(M.trace('late_mean')[:])


    

    




/projects/anaconda3/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j






    
Out[38]:





<matplotlib.axes._subplots.AxesSubplot at 0x7fbd65f84470>






    

In [39]:

    

pymc.Matplot.plot(M)


    

    




Plotting switchpoint
Plotting early_mean
Plotting late_mean






    













    













    

In [ ]: