from numpy import array
%load_ext rmagic 
First, load the data and plot it.   Note that we have years in actual years, and idx as a counter for years.
disasters = ( 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)
years = range(1851,1962)
len(disasters)
idx = range(111)

%Rpush years disasters
%R plot(years, disasters)
Next, obtain cumulative totals.   Note that this is a slightly longwinded way of doing things, but we had an ulterior motive for wanting to use as much linear algebra as possible.
import numpy as np
idxl = np.tril_indices(111)
idxu = np.triu_indices(111)
##create blank matrices of zeros
a=np.zeros(shape=(111,111))
b=np.zeros(shape=(111,111))
## Fill the upper triangle with ones
a[idxl] = 1
b[idxu] = 1
## Now get cumulative sums
cumsumafter =  np.dot(b, disasters)
cumsumbefore =  np.dot(a, disasters)
## Now append zero at start/end as appropriate.
cumsumafter = np.append(cumsumafter, 0)
cumsumbefore = np.append(0, cumsumbefore)

Very simple plot of cumulative sums seems to show a change point.
%Rpush years cumsumbefore
%R plot(years, cumsumbefore[-1])

from scipy import random, exp

#Provide values for the priors
alphaprior = 0.001
betaprior = 0.001
gammaprior = 0.001
deltaprior = 0.001





def mprob(lb, la, ap, bp, gp, dp, m, sm, sl):
        return (lb**(ap+sm-1))*(exp(-1*((bp)+m)*lb))*(la**(gp+sl-1))*(exp(-1*((dp)+(111-m))*la))
blah=np.zeros(111)
for i in range(111):
    blah[i] = mprob(lambdabefore, lambdaafter, alphaprior, betaprior, gammaprior, deltaprior, i, cumsumbefore[i], cumsumafter[i]) 
    
condalpha =  blah / sum(blah)
m=random.choice(range(111), size=1, replace = 0, p=condalpha)

print lambdabefore
print lambdaafter
outlambdabefore=np.zeros(5000)
outlambdaafter=np.zeros(5000)
outm=np.zeros(5000)

outlambdabefore[0]=12
outlambdaafter[0]=3
outm[0]=60


for iter in range(1,5000):
    outlambdabefore[iter] = random.gamma(alphaprior+cumsumbefore[outm[iter-1]], 1.0/(betaprior+outm[iter-1])) 
    outlambdaafter[iter] = random.gamma(gammaprior+cumsumafter[outm[iter-1]], 1.0/(deltaprior+111-outm[iter-1])) 

    blah=np.zeros(111)
    for i in range(111):
        blah[i] = mprob(outlambdabefore[iter], outlambdaafter[iter], alphaprior, betaprior, gammaprior, deltaprior, i, cumsumbefore[i], cumsumafter[i]) 
    
    condalpha =  blah / sum(blah)

    outm[iter]=random.choice(range(0,111), size=1, replace = 0, p=condalpha)


    

%Rpush outlambdabefore outlambdaafter outm
%R plot(c(1:5000), outlambdabefore, type = "l")
%R par(mfrow = c(2,2)); hist(outlambdabefore[c(1000:5000)], xlim = c(0,5)); hist(outlambdaafter[c(1000:5000)], xlim =c(0,5))
%R plot(xtabs(~outm[c(1000:5000)]))
%R acf(outlambdabefore[c(1000:5000)]) ; acf(outlambdaafter[c(1000:5000)])