5
$\begingroup$

Suppose I have a potential energy function $V(x)$ and I'd like to draw random samples of the position of the particle, $x$, from the Boltzmann distribution

$$P(x) = \frac 1 Z e^{-\beta V(x)}$$

at inverse temperature $\beta$. Of course, if I have a harmonic oscillator, this is just a normal distribution. But is there a method that works for a general $V(x)$?

And can the method be generalised to 2D or 3D systems - i.e., $V(x,y)$ or $V(x,y,z)$?

Improve this question
$\endgroup$
3
  • 6
    $\begingroup$ I'm not sure this is really a physics question. Some sources you might find valuable are: en.wikipedia.org/wiki/Inverse_transform_sampling for an analytic approach, en.wikipedia.org/wiki/Rejection_sampling for a straightforward numerical algorithm that works in low numbers of dimensions, en.wikipedia.org/wiki/Markov_chain_Monte_Carlo if you need more dimensions. In general, sampling from an arbitrary distribution with a lot of variables is a hard computational problem. If the distribution is of a "nice" form or has low numbers of dimensions it becomes easier to do. $\endgroup$ Commented Jan 21 at 23:34
  • 1
    $\begingroup$ I would also suggest asking this at the stats stack exchange: stats.stackexchange.com. I suspect people will point you to a previously answered question with more information than my short comment. $\endgroup$ Commented Jan 21 at 23:35
  • 1
    $\begingroup$ More generally, this is what the metropolis (or any algorithm called monte carlo something) tries to do. $\endgroup$ Commented Jan 21 at 23:43

2 Answers 2

Reset to default
10
$\begingroup$

If the normalization factor $Z$ can be easily calculated, one can use inverse transform sampling to implement the importance sampling methods, such as rejection sampling. If the normalization cannot be easily calculated, other methods must be used. The simplest is probably the rejection sampling. However, if the probability density is highly concentrated in a subset of the configuration space, the Markov Chain Monte Carlo method provides a more efficient alternative. This is the case in Statistical Mechanics, where a reliable partition function calculation is almost always hopeless due to the complex nature of the energy landscape.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ @Andrew, you're right. I wrote "rejection" method, although I had in mind the "inverse transform" method. I have corrected the text. $\endgroup$ Commented Jan 22 at 23:21
3
$\begingroup$

A simple method is to build a process whose stationary distribution is precisely your target distribution. Then sample trajectories of the process and take points of these trajectories that have enough separation to overcome correlations. These points should be distributed according to your target distribution. In more detail, for multiple dimensions, the process $$ d \vec{X}_t = -\nabla\cdot V\left(\vec{X}_t\right)\,dt+\sqrt{\frac{1}{\beta}}\,d\vec{W}_t, $$ where each component of $\vec{W}_t$ is a Wiener process (without correlations among them), has your target distribution as stationary distribution. Say $M$ is the number of samples that you want to get, then the algorithm would look something like:

  1. Initialize the process, I would select $\vec{X}_0$ as a minima of the potential. If $V$ has many different minima this could be an important step.
  2. Select integration time step $\Delta t$.
  3. Select number of iterations between measures $n$.
  4. Iterate the following formula $n\cdot M$ times: $\vec{X}_{t+\Delta t} = -\nabla\cdot V\left(\vec{X}_t\right)\,\Delta t+\sqrt{\frac{\Delta t}{\beta}}\,\vec{G}(0,1)$. Where $\vec{G}(0,1)$ is a vector of random numbers with zero mean and variance 1.
  5. Your samples correspond to the states $\{\vec{X}_{n\Delta t},\vec{X}_{2n\Delta t},\dots,\vec{X}_{Mn\Delta t}\}$.

The tricky steps are to fix $n$ and $\Delta t$ to ensure precision at reasonable computation times.

The nice thing about this method is that it is quite general, however, it is very inefficient with many dimensions (as it requires evaluating $V$ in every iteration). There are many other methods that you can use, such as rejection, inversion method, etc. The best method to use depends on your particular problem (form of $V$). This is a nice reference for methods of random number generation: Stochastic Numerical Methods: An Introduction for Students and Scientist.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.