``````

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 [ ]:

``````