# Monte Carlo sampling of probability distribution

**URL:** <https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726>\
**Category:** version agnostic\
**Tags:** modeling\
**Created:** [February 1, 2024, 1:27am UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726 "2024-02-01T01:27:13Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![m8cNzLHfUBw.FptZ](https://avatars.discourse-cdn.com/v4/letter/m/65b543/32.png) [@m8cNzLHfUBw.FptZ](https://discourse.pymc.io/u/m8cNzLHfUBw.FptZ)\
**Post date:** [February 1, 2024, 1:27am UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/1 "2024-02-01T01:27:13Z")

</div>

**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](https://discourse.pymc.io/t/sampling-from-prior-predictive-distribution/141) 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](https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/pymc_overview.html), 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!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [February 1, 2024, 8:05am UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/2 "2024-02-01T08:05:46Z")

</div>

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)

```python
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

---

<div class="post-metadata">

**Author:** ![m8cNzLHfUBw.FptZ](https://avatars.discourse-cdn.com/v4/letter/m/65b543/32.png) [@m8cNzLHfUBw.FptZ](https://discourse.pymc.io/u/m8cNzLHfUBw.FptZ)\
**Post date:** [February 1, 2024, 2:47pm UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/3 "2024-02-01T14:47:12Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [February 2, 2024, 10:01am UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/4 "2024-02-02T10:01:50Z")

</div>

> 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](https://www.pymc.io/projects/docs/en/latest/api/generated/pymc.sampling.mcmc.sample.html)

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](https://www.pymc.io/projects/docs/en/latest/api/generated/pymc.smc.sample_smc.html)

---

<div class="post-metadata">

**Author:** ![m8cNzLHfUBw.FptZ](https://avatars.discourse-cdn.com/v4/letter/m/65b543/32.png) [@m8cNzLHfUBw.FptZ](https://discourse.pymc.io/u/m8cNzLHfUBw.FptZ)\
**Post date:** [February 2, 2024, 7:02pm UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/5 "2024-02-02T19:02:30Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [February 2, 2024, 8:32pm UTC](https://discourse.pymc.io/t/monte-carlo-sampling-of-probability-distribution/13726/6 "2024-02-02T20:32:38Z")

</div>

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()])` ?