Monte Carlo sampling of probability distribution

Questions version agnostic
1

Setup

I have a Potts model with spins s_i \in \{1, ..., q\}.:

P(s_1, ..., s_n) \sim \exp(\sum_i h_i(s_i) + \sum_{i < j} J_{ij}(s_i, s_j)),

where P denotes the probability of finding the configuration (s_1, ... s_n). I am given fixed values for h_i(s_i) for each s_i and J_{ij}(s_i, s_j) for each pair (s_i, s_j) (so I am not doing Bayesian inference to learn s_i or J_{ij}).

What I would like to do

Start from a random initialization point, then run a bunch of Monte Carlo steps to equilibrate it, to get a configuration C with high probability. Then starting from that one configuration C, run many independent Monte Carlo trajectories (all starting from C) (each for a certain number of time steps).

My main question is: how can I do this in pymc?

Notes

  1. I understand that I am not using most of the machinery of pymc (i.e., Bayesian inference). I have a probability distribution that I want to sample from using Markov Chain Monte Carlo, all from the same starting point C. I want to use pymc instead of writing it myself because pymc includes optimized versions of HMC, NUTS, etc. that would take me a long time to write myself.
  2. I looked around this forum and found this post showing that I can sample from a distribution without doing the full Bayesian inference. This is great because I have no data; I just want to sample from a given Potts distribution.
  3. If I have the above working for Metropolis-Hastings, what are the steps I need to get Hamiltonian Monte Carlo to work?
  4. One question I have is how to specify a common starting point C for a bunch of independent Monte Carlo trajectories.
  5. I am new to pymc so any advice you can provide will be super helpful.
  6. Many resources I see are about using the full power of pymc (Bayesian inference), as opposed to simply doing sampling from a specified distribution. For example, many use standard distributions like Gaussian; how can I specify instead the Potts model above with given (fixed) values for h_i and J_{ij}?
  7. The reason I ask question 6 is that from the pymc overview here, it talks about how one can define new functions (which I need to since there is no pre-defined Potts models, unlike the Gaussian distribution) but it says that if one defines new functions, it is not possible to compute the gradient so I cannot use HMC or NUTS samplers. I would like to not just use Metropolis-Hastings, but I would also like to use HMC/NUTS samplers.

Thank you!

2

I don’t follow the whole starting point thing, but here is how you can use pymc to sample from a pdf (in this case a normal but could be anything)

with pm.Model() as m:
  x = pm.Flat("x")
  pm.Potential("pdf", pm.logp(pm.Normal.dist(), x))
  idata = pm.sample()  # posterior x will follow a unit normal density

I used pm.logp helper to give me the Normal density but I could have written it directly with PyTensor operations

3

Thank you! This is very helpful.

  1. If you had written the suggested solution using PyTensor operations, would one still be able to use gradient-based (i.e., HMC/NUTS) samplers? My understanding is that writing functions myself in pymc means I cannot use gradient-based samplers like HMC/NUTS.
  2. Is using pm.sample() the same as sample_prior_predictive as long as I do not include data in observed=data (in my case I do not have observed data; I just want to sample a distribution)?
  3. I’m assuming the above code will use Metropolis-Hastings or HMC; is that correct?
  4. Is there a way to specify a starting configuration where I want the sampling to start from (obviously if it is a random sampling it wouldn’t make sense, but I am using sequential Monte Carlo, so choice of starting point makes sense)? For example, let’s say I want to run the suggested code for 100 time steps so I am at a configuration C, then initialize another Metropolis-Hastings or HMC with C as a starting point. Is there a way to do that?

Thank you very much!

4

If you had written the suggested solution using PyTensor operations, would one still be able to use gradient-based (i.e., HMC/NUTS) samplers? My understanding is that writing functions myself in pymc means I cannot use gradient-based samplers like HMC/NUTS.

I wrote the suggestion with PyTensor operations. pm.logp returns an expression that uses PyTensor operations. I could have done the same manually. If you write with PyTensor operations (which should be very similar to writing numpy code) you get gradients for free and can sample with NUTS.

Is using pm.sample() the same as sample_prior_predictive as long as I do not include data in observed=data (in my case I do not have observed data; I just want to sample a distribution)?

Yes, they should be sampling from the same distribution. The distinction is prior_predictive uses forward sampling (requiring random number generators), and sample uses mcmc (requiring logp expressions and optional gradients for those)

I’m assuming the above code will use Metropolis-Hastings or HMC; is that correct?

The example code I wrote will use NUTS.

Is there a way to specify a starting configuration where I want the sampling to start from (obviously if it is a random sampling it wouldn’t make sense, but I am using sequential Monte Carlo, so choice of starting point makes sense)? For example, let’s say I want to run the suggested code for 100 time steps so I am at a configuration C, then initialize another Metropolis-Hastings or HMC with C as a starting point. Is there a way to do that?

You can define the starting point by passing initvals to sample: pymc.sampling.mcmc.sample — PyMC dev documentation

However there’s a complication to your picture. These samples need tuning before they start giving correct draws, and during tuning they will drift from your suggested initvals. If you want more fine-grained control you have to either initialize the step samplers already tuned manually, or perhaps override their tuning routine. This is non-trivial, specially for advanced samplers like NUTS.

For sequential monte carlo, you can provide the distribution of initial particle positions via the start kwarg I think, otherwise it attempts to initialize from the prior (using forward sampling): pymc.smc.sample_smc — PyMC dev documentation

5

Thank you very much!! I marked it as the solution because it answers almost everything.

I wanted to ask about a minor detail: in the suggested code above, how would I change it to use:

  1. HMC
  2. sequential Monte Carlo

Thanks!

6

SMC has a separate function, you would call pm.sample_smc instead of pm.sample.

For HMC, instead of NUTS, you have to specify the step, something like pm.sample(..., step=[pm.HMC()]) ?

1 Like