Given $3$ real variables $x_1, x_2, x_3 \equiv \bf{x}$, consider their probability density function (PDF) \begin{equation} P({\bf x}) = C \, p(x_1) \cdots p(x_3) \exp[f({\bf x})], \end{equation} where $p(x)$ is a known one-dimensional PDF, $C$ a normalization constant, and $f$ is defined as follows. Consider three binary variables, or 'spins', $S_i = \pm 1, i=1,..,N$, their Hamiltonian \begin{equation} H_{\bf x}[{\bf S}] \equiv - (x_1 S_1 S_2 + x_2 S_1 S_3 + x_3 S_2 S_3), \end{equation}
where the $x$s can be interpreted as the couplings between spin pairs. Let me denote by ${\bf S}^1_{{\bf x}}$ and ${\bf S}^{2}_{{\bf x}}$ the ground state and the first excited state of $H_{\bf x}$, respectively, and set

\begin{equation} f({\bf x}) = (H_{\bf x}[{\bf S}^1_{\bf x}]-H_{\bf x}[{\bf S}^2_{\bf x}])^2. \end{equation}

Do you know an efficient method to draw random samples from $P$, given the particular form above of $P$?

In particular, the form of $P$ above is such that random samples from the factorized distributions $p(x)$ can be easily drawn with inverse transform sampling.

I have tried the two following methods to solve my problem:

  • The reweighting method:

    1. Consider a 'temporary' random sample ${\bf x}^1_{\rm t}$, where each of the three entries of ${\bf x}^1_{\rm t}$ is drawn independently from $p$.
    2. Repeat point 1 $S\gg 1$ times and obtain $S$ temporary samples ${\bf x}^1_{\rm t}, \cdots, {\bf x}^S_{\rm t}$.
    3. Introduce the weight of each of these samples \begin{equation} w^s \equiv \frac{e^{f({\bf x}^s_{\rm t})}}{\sum_{p=1}^S e^{f({\bf x}^p_{\rm t})}} \end{equation}
    4. Reweighting: draw a random number $r$ uniformly distributed in $[0,1)$, find the value of $s$ such that \begin{equation} \sum_{p=1}^s w^p < r < \sum_{p=1}^{s+1} w^p, \end{equation} and obtain sample ${\bf X}_1 \equiv {\bf x}^s_t$.
    5. Repeat point 4 $S\gg 1$ times and obtain samples ${\bf X}^1, \cdots, {\bf X}^S$, which are distributed according to $P$.
  • the Markov Chain Monte Carlo method.

However, both methods are not efficient for my specific problem.

  • 1
    $\begingroup$ Every PDF whose support is a cartesian product $X^N$ for some $X\subseteq \mathbb R$ can be written in that form. It seems too general. Perhaps your function $f$ has some useful extra properties? $\endgroup$ Mar 12, 2017 at 0:24
  • 1
    $\begingroup$ Please see the revised question, where I wrote $f$ explicitly. $\endgroup$
    – James
    Mar 13, 2017 at 15:00
  • 1
    $\begingroup$ By the MCMC method you mean the Metropolis algorithm? $\endgroup$
    – Wuestenfux
    Mar 15, 2017 at 18:59
  • 1
    $\begingroup$ Yes. I also tried more sophisticated MCMC methods. $\endgroup$
    – James
    Mar 16, 2017 at 15:37

1 Answer 1


Your question is a bit similar to a hanging post here How can we simulate from a geometric mixture?

To cite @whuber's comment under the SE post,

Without additional assumptions, this seems unlikely. ...Suppose that associated with each $f_i$ is an interval $I_i$ on which $f_i≤1$ and $Pr_i(I_i)>1−\epsilon$, outside of which $0<f_i<\epsilon$, and $I_i\cap I_j=\emptyset$ for $i\neq j$. Then the separate generators would almost always produce values in $I_i$, but the probability of $\prod_i f_i$ could be concentrated anywhere, seemingly unrelated to the $I_i$...

But that is a very general comment when there is nothing more we know about the densities $f_i$. Now we know that $f_i(x)=p(x)$ in your case and $$H_{\boldsymbol{x}}[\boldsymbol{S}]=-(x_{1}S_{1}S_{2}+x_{2}S_{1}S_{3}+x_{3}S_{2}S_{3})=-\left(\begin{array}{ccc} x_{1} & x_{2} & x_{3}\end{array}\right)\left(\begin{array}{ccc} S_{1}S_{2}\\ & S_{1}S_{3}\\ & & S_{2}S_{3} \end{array}\right)\left(\begin{array}{c} 1\\ 1\\ 1 \end{array}\right)=-\left(\begin{array}{ccc} x_{1} & x_{2} & x_{3}\end{array}\right)\left(\begin{array}{ccc} S_{1}\\ & S_{3}\\ & & S_{2} \end{array}\right)\left(\begin{array}{ccc} S_{2}\\ & S_{1}\\ & & S_{3} \end{array}\right)\left(\begin{array}{c} 1\\ 1\\ 1 \end{array}\right)=-\left(\begin{array}{ccc} x_{1} & x_{2} & x_{3}\end{array}\right)\pi_{1}\left(\begin{array}{ccc} S_{1}\\ & S_{2}\\ & & S_{3} \end{array}\right)\pi_{2}\left(\begin{array}{ccc} S_{1}\\ & S_{2}\\ & & S_{3} \end{array}\right)\left(\begin{array}{c} 1\\ 1\\ 1 \end{array}\right)$$ where $\pi_{1},\pi_{2}$ are two permutations in $Sym(3,\mathbb{R})$ we know that all you need is to sample $(x_1,x_2,x_3)$ and $diag(S_1,S_2,S_3)$.

To sample from $(x_1,x_2,x_3)$ we may need to concern about the problem I mentioned in the beginning yet one solution after you know $P$ is to do a tranformed MCMC or dynamic MCMC. The reason why a direct MCMC failed here is probably due to $P$'s heavy-tailedness.

As how to sample a random matrix $diag(S_1,S_2,S_3)$, I think you can choose whatever appropriate random matrix sampling method.

I do not quite understand your description about reweighting. Could you explain more in details?

  • $\begingroup$ Thank you for your comment, I detailed the explanation of the reweighting method. Also, $f(x)$ is defined in terms of $H_{\bf x}[S^1_{\bf x}]$, where $S^1_{\bf x}$ depends on $x$. It follows that $\bf x$ and $S^1_{\bf x}$ are not independent, while in your comment you seem to say that they are. $\endgroup$
    – James
    Mar 16, 2017 at 16:09
  • $\begingroup$ When $x$ and $S_x$ are not independent, you need dependent sampling, search "dependent dirichlet process" for example. $\endgroup$
    – Henry.L
    Mar 16, 2017 at 16:22
  • $\begingroup$ Also, given that they are not independent, writing $H_{\bf x}[S^1_{\bf x}]$ as a product of matrices and vectors is not useful. $\endgroup$
    – James
    Mar 17, 2017 at 9:13

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.