"Injecting" values into each MCMC iteration

Questions v5
1

I want to run a weighted regression in PyMC but with a posterior distribution of weights I have precalculated on the same data.

Say I have obtained a matrix W of N observations and Q draws (i.e., draws * chains). When calculating the regression likelihood, I want to multiply the normal log-pdf of the outcome with a vector w_q (of size N) for each MCMC iteration q.

Question:
What’s the best way to pass each MCMC iteration with a new draw from that matrix
(assume Q matches and enough weights were drew, i.e. MCMC iterations \leq Q)?

Nice to have: Ideally, I would like to hack the model definition, rather than the optimizer (step function?). So it would look something like the below, where, w is some way to fool a distribution object to spit out predefined vectors from matrix W in order.

w = ...
y_obs = pm.Potential(
    "y_obs", 
    w*pm.logp(pm.Normal.dist(mu=mu, sigma=sigma), y)
)

But I’ll take any solution you may have.

Thank you

2

What I tried so far:

  • Create a CustomDist for the y_obs with a custom logp function that draws from a python generator.
  • Create a similar logp function that goes into the pm.Potential for the y_obs.
  • Create a CustomDist for w which has a generator for both logp and random.

But these seem to use only a handful (1-10) draws from the weights matrix (I keep track with an internal counter), instead of throughout the sampling.
However, my knowledge in the inner-works of pymc is very limited, I don’t always understand what I try :slight_smile:

3

You need a custom step sampler if you want to change the logp artifically at every mcmc step

4

Thanks Ricardo.

Although slightly unfortunate as the solution will probably be less generic and harder to understand :slight_smile:

I skimmed the NUTS/BaseHMC/ArrayStepShared code, it’s quite complex and there is relatively little documentation/examples on custom steps in pymc5.
Do you have a sense what the minimal change could possibly be?
Where should I intervene? on step? astep? _hamiltonian_step?
How will the custom step be incorporated in the pm.Model()?

Many thanks