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.


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

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

#Gibbs Sampling
g=function(z, c) {

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

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

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

enter image description here


Applying A Linear Transformation

The linear transformation $$ \left.\begin{array}{l} \displaystyle u=\frac12(x+y)\\ \displaystyle v=\frac{\sqrt3}2(x-y) \end{array}\right\}\iff\left\{\begin{array}{l} \displaystyle x=u+\frac{v}{\sqrt3}\\ \displaystyle y=u-\frac{v}{\sqrt3} \end{array}\right.\tag1 $$ equates $$ x^2-xy+y^2=u^2+v^2\tag2 $$ Since areas are scaled uniformly via linear transformations, so will densities.

The square $[0,1]^2$ in $(x,y)$-space is the diamond in $(u,v)$-space with vertices $$ \overbrace{\ \ \ (0,0)\ \ \ \vphantom{\frac{\sqrt3}2}}^{(0,0)},\overbrace{\left(\frac12,\frac{\sqrt3}2\right)}^{(1,0)},\overbrace{\ \ \ (1,0)\ \ \ \vphantom{\frac{\sqrt3}2}}^{(1,1)},\overbrace{\left(\frac12,-\frac{\sqrt3}2\right)}^{(0,1)}\tag3 $$

Using A Radial Density

In $(u,v)$-space, the density is a simple radial density and we can use tools similar to the one-dimensional case, then map back to $(x,y)$-space.

In an wedge with aperture $\theta$, the cumulative weight for $u^2+v^2$ is $$ \int_0^r\overbrace{t^2}^{\text{density}}\overbrace{\theta t}^{\substack{\text{arc}\\\text{length}}}\,\overbrace{\mathrm{d}t}^{\substack{\text{radial}\\\text{width}}}=\frac\theta4r^4\tag4 $$ Thus, the cumulative density in a wedge grows proportional to $r^4=\left(u^2+v^2\right)^2$.

Generate And Discard Points

If we don't mind generating more points than we are going to use, we can generate $(u,v)$ points in a third of a circle of radius $1$ and keep only the points mapped into $[0,1]^2$ in $(x,y)$-space, which is the diamond listed in $(3)$:

enter image description here

Although the lighter gray region is much smaller than the diamond containing the points to keep, only $\frac{5\sqrt3}{4\pi}=68.916\%$ of the points generated are kept because the density is less near the origin.

With $(s,t)$ picked uniformly in $[0,1]^2$, generate $$ \bbox[5px,border:2px solid #CA0]{(x,y)=\overbrace{t^{1/4}\left(\cos\left(\frac{2\pi}3s-\frac\pi3\right),\sin\left(\frac{2\pi}3s-\frac\pi3\right)\right)}^{\large(u,v)}\begin{bmatrix}1&1\\\frac1{\sqrt3}&-\frac1{\sqrt3}\end{bmatrix}}\tag5 $$ and discard any points with $x\gt1$ or $y\gt1$.

Here is a plot of $100,\!000$ points generated with $(5)$ where only $68,\!830$ were kept (approximately $68.916\%$). The red contours are where $x^2-xy+y^2\in\left\{\frac1{64},\frac1{16},\frac9{64},\frac14,\frac{25}{64},\frac9{16},\frac{49}{64}\right\}$

enter image description here

Generate Only The Points Needed

If we don't want to discard almost $1/3$ of the points we generate, we can generate only points in the diamond, but the formulas will be more complicated than $(5)$. It is simpler to work with the upper half of the diamond (which gives the lower right half of $[0,1]^2$ in $(x,y)$-space). For the lower half (which gives the upper left half of $[0,1]^2$ in $(x,y)$-space), we can use the same formula and simply negate $v$.

In the diagram below, the polar equation for the non-radial line bounding the upper half of the diamond is $$ r=\frac{\sqrt3}2\sec\left(\theta-\frac\pi6\right)\tag6 $$ where $0\le\theta\le\frac\pi3$.

enter image description here

The cumulative weight as a function of $\theta$ is proportional to the integral of $r^4$: $$ \frac{9\sqrt3}{20}\int_0^\theta\sec^4\left(\phi-\frac\pi6\right)\,\mathrm{d}\phi=\frac{9\sqrt3}{20}\left(\tan\left(\theta-\frac\pi6\right)+\frac13\tan^3\left(\theta-\frac\pi6\right)+\frac{10\sqrt3}{27}\right)\tag7 $$ The constant $\frac{9\sqrt3}{20}$ was chosen so that when $\theta=\frac\pi3$, the total weight is $1$.

To choose a random $\theta$, we need to invert $(7)$. To do this, we make use of the formula $$ \sinh(3x)=3\sinh(x)+4\sinh^3(x)\tag8 $$ which, when multiplied by $\frac23$, becomes $$ \frac23\sinh(3x)=2\sinh(x)+\frac13(2\sinh(x))^3\tag9 $$ Applying $(9)$ to $(7)$ gives $$ \bbox[5px,border:2px solid #CA0]{\theta=\frac\pi6+\tan^{-1}\left(2\sinh\left(\frac13\sinh^{-1}\left(\frac{10}{3\sqrt3}s-\frac{5\sqrt3}9\right)\right)\right)}\tag{10} $$ where $s$ is picked uniformly in $[0,1]$.

Now that we have $\theta$, we can apply $(4)$ and $(6)$ to get $$ \bbox[5px,border:2px solid #CA0]{(x,y)=\overbrace{\left(\frac{9|2t-1|}{16}\right)^{\!\!1/4}\!\sec\left(\theta-\frac\pi6\right)(\cos(\theta),\operatorname{sgn}(2t-1)\sin(\theta))}^{\large(u,v)}\begin{bmatrix}1&1\\\frac1{\sqrt3}&-\frac1{\sqrt3}\end{bmatrix}}\tag{11} $$ where $t$ is picked uniformly in $[0,1]$.

Here is a plot of $68,\!916$ points generated with $(10)$ and $(11)$ where all points were used. The red contours are where $x^2+xy+y^2\in\left\{\frac1{64},\frac1{16},\frac9{64},\frac14,\frac{25}{64},\frac9{16},\frac{49}{64}\right\}$

enter image description here

  • $\begingroup$ hmm I like this quite a lot as it explains how to use transformations. However since my goal is to have function that randomizes data based on an input function this moves the problem for me from randomizing to "how to linearize in a generic sense", which seems equally difficult to generalize? $\endgroup$ – paul23 yesterday
  • $\begingroup$ After the simplification from an elliptically-symmetric distribution to a circularly-symmetric distribution, this answer shows two approaches to generating a two dimensional distribution along the same lines as you mentioned for the one dimensional distribution: compute the cumulative distribution and invert. Both approaches here do this by inverting the cumulative distribution for a given angle, then inverting the cumulative radial distribution. This could be done with any coordinates, polar, rectangular, etc. $\endgroup$ – robjohn yesterday
  • $\begingroup$ If the question was not how to simulate the given distribution from a uniform distribution on $[0,1]^2$, then I misunderstood and I apologize. $\endgroup$ – robjohn yesterday
  • $\begingroup$ the question states that, but I not realise it's a bit of an xy problem. The function was only given as example function that is "not trivially linearized". I realise now that this might have been doing from my part. But the goal of the question was to find a way where I have any normalized function/surface. (I hope that by looking at a specific solution I could extract a generic one got any function) $\endgroup$ – paul23 1 hour ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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