# Sampling from a particular multivariate probability distribution

Given $3$ real variables $x_1, x_2, x_3 \equiv \bf{x}$, consider their probability density function (PDF) $$P({\bf x}) = C \, p(x_1) \cdots p(x_3) \exp[f({\bf x})],$$ 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 $$H_{\bf x}[{\bf S}] \equiv - (x_1 S_1 S_2 + x_2 S_1 S_3 + x_3 S_2 S_3),$$
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

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

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 $$w^s \equiv \frac{e^{f({\bf x}^s_{\rm t})}}{\sum_{p=1}^S e^{f({\bf x}^p_{\rm t})}}$$
4. Reweighting: draw a random number $r$ uniformly distributed in $[0,1)$, find the value of $s$ such that $$\sum_{p=1}^s w^p < r < \sum_{p=1}^{s+1} w^p,$$ 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.

• 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? Mar 12 '17 at 0:24
• Please see the revised question, where I wrote $f$ explicitly. Mar 13 '17 at 15:00
• By the MCMC method you mean the Metropolis algorithm? Mar 15 '17 at 18:59
• Yes. I also tried more sophisticated MCMC methods. Mar 16 '17 at 15:37

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?

• 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. Mar 16 '17 at 16:09
• When $x$ and $S_x$ are not independent, you need dependent sampling, search "dependent dirichlet process" for example. Mar 16 '17 at 16:22
• 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. Mar 17 '17 at 9:13