Search code examples
pythonmontecarlopymc

How to monitor the proposal covariance in PyMC?

I'm testing the AdaptativeMetropolis step method in PyMC (as documented here), and wish to see it in action. Such a step method consists in block-updating some variables using a multivariate jump distribution whose covariance is tuned during sampling.

It is possible to print the proposal covariance matrix vs time? It looks to me that this is not documented.

Thanks.

EDIT: as a toy model, let us consider the case of two similar dice (let's say that they come from the same factory). We wish to estimate whether they are biased. 

import pymc
n = 100
alpha = pymc.Gamma('alpha', alpha=1, beta=1)
P1 = pymc.Beta('p1', alpha=alpha, beta=2)
P2 = pymc.Beta('p2', alpha=alpha, beta=2)
Y1 = pymc.Binomial('y1', n=n, p=P1, value=70, observed=True)
Y2 = pymc.Binomial('y2', n=n, p=P2, value=50, observed=True)
model = pymc.Model([alpha, P1, P2, Y1, Y2])
MC = pymc.MCMC(model)
MC.use_step_method(pymc.AdaptiveMetropolis, [P1, P2])

I now expect/guess that P1 and P2 are set to be sampled using a normal bivariate proposal with adaptative covariace. How can I monitor this? The sampling is done after calling

MC.sample(iter=11000, burn=10000)
Solution

  • Digging the source code, I found that it is possible to tune the verbosity to print the covariance matrix:

    if self.verbose > 1:
  print_("\tUpdating covariance ...\n", self.C)
  print_("\tUpdating mean ... ", self.chain_mean)

    Hence, it should be sufficient:

    MC.use_step_method(pymc.AdaptiveMetropolis, [P1, P2], verbose=2)

    Also the final correlation matrix is stored into MC.get_state()['step_methods']['AdaptiveMetropolis_p1_p2']['proposal_sd']. This can be used to compute the covariance matrix after the sampling (np.dot(corr, corr.T)).