How to transform a multi (2) dimensional uniform random function to a given probability density function? Well as per title, say I have the probability density function on domain $x \in [0,1] ; y \in [0,1]$
$$f(x,y) = \frac{12}{5} \left( x^2 + y^2 - xy \right)$$
Can I generate this density function from a given uniform (pseudo) random function on the same domain?
When using a single variant it's slightly easy:

*

*Integrate the function to calculate the cumulative distribution function

*calculate the inverse of the CDF.

*plug in the uniform random function.

However in multiple dimensions this can't be really done the "inverse" isn't clearly defined. - If I could split the variables it's a bit more trivial. But how can this be done in the generic case where the variables aren't independent?
I could of course do it by rasterizing the function and getting linearizing the raster (just putting row behind row) and then using normal technologies for this. However this numerical approach seems inexact and arbitrary.
 A: Yes, use the Metropolis-Hastings algorithm with a uniform proposal distribution. Gibbs sampling is a special case of Metropolis-Hastings that can also be used. They fall under the category of Markov Chain Monte Carlo and can be used to simulate from high dimensional distributions (including 2 here).
Here is working code that uses (1) Gibbs sampling and (2) Metropolis-Hastings via an independence sampler.
The Gibbs sampler works by starting by fixing $y$. It then generates from the distribution that is $f(x,y)$ with $y$ plugged in, which multiplied by a normalizing factor gives the conditional pdf $g(x|y)$. Then, the $x$ value received from this is used to generate a $y$ from $g(y|x)$. Acceptance-rejection sampling  is used to generate from the conditional pdf's. This is repeated 10,000 times and the burn-in of 100 is discarded.
The MH algorithm here uses an independence sampler, wherein the new state of the Markov chain does not depend on the previous state and is simulated uniformly from $[0,1]\times [0,1]$ each iteration.
f=function(x, y) {
  return(12/5*(x^2+y^2-x*y))
}

sample=data.frame(matrix(0, nrow=0, ncol=2))
x=c(1/2, 1/2)
for (i in 1:10000) {
  xpx=runif(1)
  xpy=runif(1)
  a=f(xpx, xpy) / f(x[1], x[2])
  r=min(1, a)
  u=runif(1)
  if (u<r) {
    x=c(xpx, xpy)
  }
  sample[i,]=x
  print(i)
}
sample_no_burnin=sample[100:10000,]
library(ggplot2)

#Gibbs Sampling
g=function(z, c) {
  return(12/5*(z^2+c^2-z*c)*5/(12*c^2+6*c+4))
}

samp=function(c) {
  while (1>0) {
    x=runif(1)
    u=runif(1)
    if (g(x, c) / 1 >= 13/5*u) {
      return(x)
    }
  }
}

y=1/2
my_sample=data.frame(matrix(0, nrow=0, ncol=2))
for (i in 1:10000) {
  x=samp(y)
  y=samp(x)
  my_sample[i,]=c(x, y)
  print(i)
}
my_samp_no_burnin=my_sample[100:10000,]

library(dplyr)
samps_no_burnin=bind_rows(my_samp_no_burnin, sample_no_burnin) %>% 
  mutate(method=rep(c("gibbs", "metropolis-hastings"), each=9901))
ggplot(samps_no_burnin) + stat_density_2d(aes(X1, X2)) + theme_bw() + facet_wrap(vars(method))

Does either look better? It's hard to say...

