fonnesbeck
3/25/2010 - 12:07 AM

## Simple Bayesian hierarchical model in PyMC. Coded by A. MacNeil.

Simple Bayesian hierarchical model in PyMC. Coded by A. MacNeil.

"""
Simple hierarchical linear model of price as a function of age.
"""

from pymc import *
from numpy import array

# Data
age = array([13, 14, 14,12, 9, 15, 10, 14, 9, 14, 13, 12, 9, 10, 15, 11, 15, 11, 7, 13, 13, 10, 9, 6, 11, 15, 13, 10, 9, 9, 15, 14, 14, 10, 14, 11, 13, 14, 10])
price = array([2950, 2300, 3900, 2800, 5000, 2999, 3950, 2995, 4500, 2800, 1990, 3500, 5100, 3900, 2900, 4950, 2000, 3400, 8999, 4000, 2950, 3250, 3950, 4600, 4500, 1600, 3900, 4200, 6500, 3500, 2999, 2600, 3250, 2500, 2400, 3990, 4600, 450,4700])/1000.

"""Original WinBUGS model:

model
{
for( i in 1 : N ) {
mu[i] <- a+b*age[i]
price[i] ~ dnorm(mu[i], tau)
}
tau ~ dgamma(0.001, 0.001)
sigma <- 1 / sqrt(tau)
a ~ dnorm(0, 1.0E-12)
b ~ dnorm(0, 1.0E-12)
}"""

# Constant priors for parameters
a = Normal('a', 0, 0.0001)
b = Normal('b', 0, 0.0001)

# Precision of normal distribution of prices
tau = Gamma('tau', alpha=0.1, beta=0.1)

@deterministic
def mu(x=age, a=a, b=b):
# Linear age-price model
return a + b*x

# Sampling distribution of prices
p = Normal('p', mu, tau, value=price, observed=True)