Simplest linear model with PyMC

Question

Say I try to estimate the slope of a simple y= m * x problem using the following data:

x_data = np.array([0,1,2,3])
y_data = np.array([0,1,2,3])

Clearly the slope is 1. However, when I run this in PyMC I get 10

slope  = pm.Uniform('slope',  lower=0, upper=20)

@pm.deterministic
def y_gen(value=y_data, x=x_data, slope=slope, observed=True):
  return slope * x

model = pm.Model([slope])
mcmc  = pm.MCMC(model)
mcmc.sample(100000, 5000)

# This returns 10
final_guess = mcmc.trace('slope')[:].mean()

but it should be 1!

Note: The above is with PyMC2.

Answer 1

You need to define a likelihood, try this:

import pymc as pm
import numpy as np

x_data = np.linspace(0,1,100)
y_data = np.linspace(0,1,100)

slope  = pm.Normal('slope',  mu=0, tau=10**-2)
tau    = pm.Uniform('tau', lower=0, upper=20)

@pm.deterministic
def y_gen(x=x_data, slope=slope):
  return slope * x

like = pm.Normal('likelihood', mu=y_gen, tau=tau, observed=True, value=y_data)

model = pm.Model([slope, y_gen, like, tau])
mcmc  = pm.MCMC(model)
mcmc.sample(100000, 5000)

# This returns 10
final_guess = mcmc.trace('slope')[:].mean()

It returns 10 because you're just sampling from your uniform prior and 10 is the expected value of that.

Answer 2

You need to set value=y_data, observed=True for the likelihood. Also, a minor point, you don't need to instantiate a Model object. Just pass your nodes (or a call to locals()) to MCMC.