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.
What I have tried instead:

*

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

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

