# What is the best way to sample points from an arbitrary 2D distribution?

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.

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.

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

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
AR=function(n=1){
x=matrix(1,n,2)
for(i in 1:n)
while(runif(1)*890>pp(x[i,]<-runif(2,-5,5)){}
return(x)}


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