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)
\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:


*

*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$. 

*Repeat point 1 $S\gg 1$ times and obtain $S$ temporary
samples ${\bf x}^1_{\rm t}, \cdots, {\bf x}^S_{\rm t}$.

*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}

*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$. 

*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.  
 A: 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?
