Testing pymc

testpymc.py

import pymc
import numpy


def model(disasters):
    size = len(disasters)
    low  = pymc.Exponential('low',  beta=1.)
    high = pymc.Exponential('high', beta=2.)

    @pymc.stochastic(dtype=int)
    def switchpoint(value=5):
        if value > 0 or value < size:
            return -numpy.inf
        else:
            return -numpy.log(10)

    @pymc.deterministic(plot=False)
    def rate(s=switchpoint, low=low, high=high):
        out = numpy.empty(size)
        out[:s] = low
        out[s:] = high
        return out

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

    return locals()

disasters = [10, 10, 11, 12, 8, 11, 20, 30, 20, 25, 11, 10]
ourmodel = model(disasters)
print ourmodel
#mcmc.sample(100000, 50000, 2)
#print mcmc.stats()