BayesFactor

#!/usr/bin/env python
# encoding: utf-8
"""
bayesfactor.py

Computation of Bayes factor for comparing two models, following the example in 
panel 7.1 in Link and Barker (2010). The posterior model probabilities are the 
product of the BF and the prior probabilities, so the BF can be calculated by 
dividing the posterior odds (the inverse of the complement of the mean of 
'true_model' int the model below) by the prior odds.

Created by Christopher J. Fonnesbeck on 2010-12-15.
Copyright (c) 2010 Vanderbilt University. All rights reserved.
"""

from pymc import Uniform, Lambda, observed, Bernoulli, geometric_like, poisson_like

# Data
Y = [0,1,2,3,8]

# Prior model probabilities
pi = (0.1, 0.9)

# Index to true model
true_model = Bernoulli('true_model', p=pi[1], value=0)

# Poisson mean
mu = Uniform('mu', lower=0, upper=1000, value=4)

# Geometric probability
p = Lambda('p', lambda mu=mu: 1/(1+mu))

@observed
def Ylike(value=Y, mu=mu, p=p, M=true_model):
    """Either Poisson or geometric, depending on M"""
    if M:
        return poisson_like(value, mu)
    return geometric_like(value+1, p)