aloctavodia/doing_bayesian_data_analysis

Python/PyMC3 versions of the programs described in Doing bayesian data analysis by John K. Kruschke

Doing_bayesian_data_analysis

This repository contains the Python version of the R programs described in the great book Doing bayesian data analysis (first edition) by John K. Kruschke (AKA the puppy book).

All the code is adapted from the Kruschke's book, except hpd.py that is taken (without modifications) from the PyMC project.

The name of the programs are the same used in the book, except they begin with a number indicating the chapter. All programs are written in Python and instead of BUGS/JAGS the PyMC3 module is used.

Thanks to Brian Naughton the code is also available as an IPython notebook

Second edition

If you are interested on the PyMC3 code for the second edition of Doing bayesian data analysis, please check this Repository.

doing_bayesian_data_analysis's Issues

Two way anova

Hi! I was wondering if you managed to get the Salary program (two way anova) working. I'm trying really hard but I don't seem to get an answer similar to the book. (for example there's no main effect on FT1 FT2 and FT3, plus so many convergence problems)

make_indices is not definied

In hpd.py on line 49 you have function make_indices which you did not defined, or imported.

09_FilconPyMC.py -- Choice of step method

Why is NUTS used as the step method for kappa? Using Metropolis appears to produce the same results, but runs faster.

I am new to pymc3 and trying to understand when to use NUTS.

Linear regression (chapter 16)

https://github.com/aloctavodia/Doing_bayesian_data_analysis/blob/master/16_SimpleLinearRegressionPyMC.py

Hi! I've tried to do linear regression by myself, using the following model (very much similar to the book), but no matter what I do, I can't manage to get a straight fit.

I've seen that in your code you do a transformation. Why is that? Is there any way to get the result without resorting to the transformation?

import numpy as np
import pymc as pm

true_intercept = 10
true_slope = 0.5 

x = np.linspace(0, 10, 100)
y = true_slope * x + true_intercept
# Jitter x and y
x += 0.5 * (np.random.rand(len(x)) * 2.0 - 1.0)
y += 1.0 * (np.random.rand(len(x)) * 2.0 - 1.0)
scatter(x, y)
show()
# Linear regression
b0 = pm.Normal('b0', 0.0, tau=1e-5, value=0)
b1 = pm.Normal('b1', 0.0, tau=1e-5, value=0)

@pm.deterministic
def mu(value=x, b1=b1, b0=b0):
    return x * b1 + b0

tau = pm.Gamma('tau', 0.01, 0.01)

resp = pm.Normal('y', mu, tau, value=data, observed=True)

mcmc = pm.MCMC([b0, b1, mu, tau, resp])
mcmc.sample(40000, burn=20000, thin=10)