I want to sample points $(x,y)$ randomly according to the Himmelblau function

$$f(x,y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2\qquad -5\le x,y\le 5$$

which I treat as a multivariate probability density function. A visualization of the function can be found here. To put it simply, what I need in the end is a collection of points which are distributed such that they resemble the Himmelblau function in the sense that I end up with more points lying near the maxima than the minima of the function.

However I am struggling to find a good way to accomplish this in practice. As far as I understand inverse transform sampling is not applicable here.

What I have tried instead:

  1. Rejection method: The function has the range of values (-5,5), (-5,5), ~800), so I tried to sample 3D points $(x,y,z)$ uniformly in this space and if $z < f(x,y)$ I would keep $(x,y)$ as my sampled point. However, this method was computationally expensive and therefore infeasible.
  2. Using numpy's choice function: I created a discrete grid of values (-5,5) x (-5,5) and for every point $(x,y)$ in this grid I computed $f(x,y)$ and stored it as the value of the grid. Then I flattened the grid to a long 1D array and used numpy's choice function to sample points from this array by using the function values as weights. For every of the sampled points I then recomputed the $(x,y)$ values. This seemed to work quite well, but I am not sure if this approach really leads to correct results.

So I was wondering if there exists another, more straightforward method for sampling points from an arbitrary distribution. Thank you.

  • 3
    $\begingroup$ 1. I do not understand why the rejection method is infeasible. $\endgroup$
    – Xi'an
    Feb 20, 2021 at 16:00
  • $\begingroup$ 2. The cdf inversion method only works in dimension one. $\endgroup$
    – Xi'an
    Feb 20, 2021 at 16:00
  • 1
    $\begingroup$ 3. A standard for simulating $(X,Y)$ in general settings is to simulate $X$ from the marginal then $Y|X$ from the conditional. $\endgroup$
    – Xi'an
    Feb 20, 2021 at 16:01

1 Answer 1


Here are two workable solutions for this specific density, rather than for an "arbitrary distribution" as requested by the title.

1. Since (thanks to Wolfram integrator!) \begin{align}\int_{-5}^z \{(x^2+y−11)^2&+(x+y^2−7)^2\}\,\text dx=(5 + z) (\\ &15 y^4 + 15 y^2 (-18 + z) + 10 y (-8 - 5 z + z^2)\\ & + 3 (775 + 15 z - 10 z^2 - 5 z^3 + z^4)))/15\\ &= G(z|y)\end{align} $$\int_{-5}^5 \{(x^2+y−11)^2+(x+y^2−7)^2\}\,\text dx= 10 (360 - 16 y - 39 y^2 + 3 y^4)/3=\mathfrak g(y)$$ and $$F(z)=\int_{-5}^z 10 (360 - 16 y - 39 y^2 + 3 y^4)/3\,\text dy\\ \qquad = 2 (11250 + 1800 z - 40 z^2 - 65 z^3 + 3 z^5)/3$$ the normalising constant of the density is$$\mathfrak z= 3/41000$$ Simulating $Y$ can thus be conducted by inverting numerically$$F(Y)=\mathfrak z U\sim\mathcal U(0,1)$$and given a realisation $y$ inverting (numerically)$$G(X|y)=\mathfrak g(y)U\sim\mathcal U(0,1)$$

2. However, a direct Accept-Reject strategy works as well. Since $$\max_{x,y} p(x,y) = \max_{x,y} (x^2+y−11)^2+(x+y^2−7)^2 = 890$$ simulating $(X,Y)\sim\mathcal U_{(-5,5)^2}$, $U\sim\mathcal U_{(0,890)}$, and accepting the simulation when $U\le p(X,Y)$ returns simulations from $p$. Here is an illustration in R:

pp=function(x)(-11 + x[1]^2 + x[2])^2 + (-7 + x[1] + x[2]^2)^2
  for(i in 1:n)

as shown by the picture below, with 500 realisations plotted on top of the log-density (Warning, the low and high regions of the target are relative, meaning there is still a considerable amount of mass near the lowest points):

enter image description here


Your Answer

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

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