Search code examples
correlationbayesianpymc

Bayesian Correlation using PyMC

I'm pretty new to PyMC, and I'm trying to implement a fairly simple Bayesian correlation model, as defined in chapter 5 of "Bayesian Cognitive Modeling: A Practical Course", which is as defined below:

Graphical model

I've put my code in an ipython notebook here, a code snippet is as follows:

mu1 = Normal('mu1', 0, 0.001)
mu2 = Normal('mu2', 0, 0.001)
lambda1 = Gamma('lambda1', 0.001, 0.001)
lambda2 = Gamma('lambda2', 0.001, 0.001)
rho = Uniform('r', -1, 1)

@pymc.deterministic
def mean(mu1=mu1, mu2=mu2):
    return np.array([mu1, mu2])

@pymc.deterministic
def precision(lambda1=lambda1, lambda2=lambda2, rho=rho):
    sigma1 = 1 / sqrt(lambda1)
    sigma2 = 1 / sqrt(lambda2)
    ss1 = sigma1 * sigma2
    ss2 = sigma2 * sigma2
    rss = rho * sigma1 * sigma2
    return np.power(np.mat([[ss1, rss], [rss, ss2]]), -1)

xy = MvNormal('xy', mu=mean, tau=precision, value=data, observed=True)

M = pymc.MCMC(locals())
M.sample(10000, 5000)

The error I get is "error: failed in converting 3rd argument `tau' of flib.prec_mvnorm to C/Fortran array"

I only found one other reference to this error (in this question) but I couldn't see how to apply the answer from there to my code.

Solution

  • This uninformative error is due to the way you have organized your data vector. It is 2 rows by n columns, and PyMC expects it to be n rows by 2 columns. The following modification makes this code (almost) work for me:

    xy = MvNormal('xy', mu=mean, tau=precision, value=data.T, observed=True)

    I say almost because I also changed your precision matrix to not have the matrix power part. I think that the MvNormal in your figure has a variance-covariance matrix as the second parameter, while MvNormal in PyMC expects a precision matrix (equal to the inverse of C).

    Here is a notebook which has no more errors, but now has a warning that requires additional investigation.