2

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.

6
  • I have no idea neither about Python nor Bayesian statistics, but this seems like a runtime problem somewhere deep in the code. You can probably do nothing about it, but you may want to report the issue to the package developer(s) and seek their advice.
    – StasK
    Dec 16, 2014 at 14:27
  • The problem with migrating it, is that after Abraham D Flaxman's partial answer, this is now definitely a Bayesian modelling issue, and not a coding issue ... unlikely to get an answer here.
    – tdc
    Dec 17, 2014 at 12:30
  • 1
    Argh. The developers should have provided an informative error message. That an input matrix has wrong dimensions is something very easy to check, and they should have done this.
    – StasK
    Dec 18, 2014 at 16:03
  • I guess this might not be relevant anylonger, but I have found your main issue. Instead of ss1 = sigma1 * sigma2 you should write ss1 = sigma1 * sigma1. Unfortunately, the results still are not "correct" all the time somehow.
    – fsociety
    Aug 19, 2015 at 12:21
  • @barsch good spot! I've updated the notebook to reflect this. The rho values seem to be pretty wild as you say.
    – tdc
    Aug 19, 2015 at 14:44

1 Answer 1

3

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.

4
  • I've marked it as accepted, since you fixed the bugs (thanks!), but I'm still not getting reasonable results
    – tdc
    Dec 18, 2014 at 13:45
  • try another question, then... this should be do-able. To address the warnings, I would look at the initial values for lambda. Also ss1 might not be defined as you intended, and the parameter for mu is precision, not variance. Dec 19, 2014 at 15:43
  • Thanks @AbrahamDFlaxman The link to the notebook doesn't work or stopped working. Any chance we can bring it up again? Dec 24, 2014 at 13:40
  • maybe a glitch in nbvviewer? it is working for me now. Dec 24, 2014 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.