Uniform Distributions Part 2: Beware of Transformations

Ok, now getting to the post that I actually meant to do to start the series on uniformity. The previous post was more about the thrilling adventure of bringing together a bunch of disparate math tools to resolve a seemingly simple question about uniform random variables. It was, ironically, not very visual, but it still in the spirit of the blog, to encourage the attitude that experiment, exploration, and trying different approaches to things is as much a part of math as it is for other scientific endeavors. 

Today we bring something much more visual and play around a lot with random number generation. It is a modernization and generalization of this post (it was the featured post on NestedTori for a looong time). In some sense, this might have more fun submission SoME, but given the time constraints (and the admitted surprise of it being earlier this year!), I opted for the adventure story that I had a much more solid narrative for, and something much more self-contained (and in fact, a requirement for the competition). We'll see how I do, but it speaks to the kind of choices and gambles one takes when choosing something to present to a larger event.

But anyway, in this post, we work toward obtaining uniform distributions of spheres and balls (this post will still be primarily setup on the disk in the plane; I'm splitting the already written post up in an effort to cap the length of my posts). We're working up to exploring a somewhat surprising way that uniformity can propagate itself. There has been some online buzz about balls and spheres and their very interesting volume formulas involving one of our frequent friends, the Gamma function, our factorial generalizer. Actually, I'd wanted to do a post about that formula for years, ever since the Pi vs. Tau debate (more specifically: volumes of spheres and balls are where it seems $\pi$ is better), I'm glad it got some visual stuff eventually, but of course, if I wanted to do my own post, maybe it would have more impact if I could rework it in some way in line with the theme of this blog. And that way, as we have been really enjoying lately, is through the lens of probability. However, we will not yet get to the full generalization today, as there is considerable material to get through (considerable FUN material, of course, but still, yours truly does admittedly have a bit of a problem keeping things short).

Summing the Squares: Not So Triangular?

First off, if we ask ourselves the question, what is the distribution of the sum of the squares of two uniform random variables, which is different from just adding them together to form the triangular distribution. This tells us what is the squared distance from the origin in the plane. Two uniform variables in $[-1, 1]$ mean that this radial squared function actually has a very complicated dependence because it has a different value in the corners than on the sides. We'd expect it to be uniform up to the unit disk, touching the boundary, and then it behaves differently when it goes into the corners. To see the real deal, we need to compute the distribution of squares of a uniform $[-1, 1]$ variables, and adding two independent instances of these. This, just as in the un-squared case, uses convolution.

So first let's examine the square of a uniform variable (first on $[0,1]$; the case $[-1,1]$ is a rescaling but causes some complications with having two signs in a square root): the pdf is $dx$ (keeping pdfs multiplied by $dx$ helps me remember how it transforms. One read/understand them more formally as infinitesimals—a "little bit of x"—or differential forms, my favorite mathematical tool in the known universe), but then taking $y = x^2$ means $x = \sqrt{y}$ and $$dx = \frac{1}{2\sqrt{y}} dy,$$ on the interval $[0,1]$ (and zero elsewhere). Write the resulting random variable $Z$ and $\rho_Z(z) = \frac{1}{2\sqrt{z}}$. To get the sum of two independent versions (call it $S$) of this requires us to convolve this function with itself. This turns out to be one tricky problem. First, sums of independent random variables are given by convolution (a fact we have used before):

\[\rho_{S}(\xi) = \int_{\mathbb{R}} \rho_{Z}(\xi - z) \rho_{Z}(z)\; dz\]

Looks straightforward enough to plug in the formula, but the trickiness comes in determining the right bounds of integration. Central to it is to reflect one of them and determine three different stages of overlap. There's 4 stages: no overlap with the reflected-translated function way to the left (so no overlap and so the function is $0$), overlap of the two unbounded sections, overlap of the two bounded sections, and no overlap due to the reflected-translated being far to the right. It's not the most inspiring calculation, so let me just show the result:

\[\rho_Y(y) = \frac{1}{2} \sin^{-1}\left(\frac{2}{y} -1\right)\]

pdf of the sum of two squared uniform random variables

Testing it with random number generators (see accompanying notebook), with a million values of $X$ and $Y$, we similarly bin the squared distance values $S$:

Very cool! You never would have thought of that shape arising from such a simple process. However, if you think about it a bit, this directly reflects that the disk is uniformly covered, but going into the corners is more difficult. The sharp break occurs because the area of the circle hits the boundary square and with distance > 1 cannot be achieved in the middle of the sides, but can be achieved at the corners. 

Now what we do is get rid of the samples $> 1$ which is basically tossing out the fraction $1 - \frac{\pi}{4}$ of the ones that exceed 1. The expected number of trials we need to run to get that is $4/\pi$ (as it's a geometric random variable), so we expect to have $\frac{\pi}{4} \approx 0.7853$ proportion of the generated samples to satisfy that condition. Then now we get to the good part: we can restrict to where $S \leq 1$. Let's just rename it to be $S$ even though technically it is different variable, with a restriction. In our random sample of a million, it did indeed leave us 784560 points (in the notebook version of this, you may get something slightly different; all that is important is it is close to $1-\pi/4$ of a million).

Now if we go back to our random variables $X$ and $Y$ (going from $[-1, 1]$ but whose squares, summed to $S$, gives a distribution exactly like the above by symmetry), also now both restricted so as to never be greater distance than $1$, we scatterplot the random samples corresponding to what produced the above:

Does it look uniform?

A slight lie: we actually showed only 1/100 of the sampled points (about 7845 in this image) because a million is wayy dense for an image like this. Looks very uniform! Again, we should emphasize the fact that although the original random variables for the $x$ and $y$ coordinates were independent uniform random variables their restrictions to the disk, $X$ and $Y$, are NOT independent. It is intuitive to understand at a big-picture, global level: $X$ and $Y$ cannot be independent, because the value of $Y$ is constrained to a very specific range once a value of $X$ is picked. The reason why this can be confusing is because when you look at an integral, the integrand $dx\; dy$ looks like a product. If we explicitly write an indicator function $\chi_D(x, y)$ then it becomes more clear, but often a double integral is written with the domain as a subscript, which hides the fact that there is a nontrivial dependence between the variables. We'll look at this in detail shortly.

Beware of Transformations

Of course, when we routinely dealing with something with circular symmetry, like a disk, often the concept of polar coordinates comes to mind. It would seem intuitive that every radius from $0$ to $1$ is equally possible, and every angle in $[0, 2\pi]$ is, also. But let's take a look at what happens:

Oops! All that uniformity in $\theta$ is scrunched down near the origin, making the sampling very concentrated here. The intuition is that the polar coordinate transformation takes what is evenly spread out in the rectangle $[0,1] \times [0,2\pi]$ (imagine it as a grid of rectangles of equal area), and then warps a bunch of stuff, so that the transformed rectangles have different sizes. This demonstrates exactly why we need Jacobians when we deal with transforming the domains of integration.  To remedy this situation, recall that the area element of the disk is not $\frac{1}{2\pi} dr \; d\theta$ which is effectively what we have here, but rather $\frac{1}{\pi}r \; dr \; d\theta$ instead. The solution to understand things in these terms is to factor it $\left(2r\; dr \right)\left(\frac{1}{2\pi} d\theta\right)$. The CDF for $r$ is simply $r^2$, so we use its inverse n the generation of $r$ from a uniform random variable. A variant of this fix is essentially what we did in our previous post on uniform distribution on the sphere. The moral of the story is, just because we start out with uniform (independed) variables, it doesn't mean we can apply an arbitrary transformation to them and still get something uniform. We may have to tweak one or more of the variables.

Breaking Down Multidimensional Uniformity via Conditioning

We explore another way of generating multidimensional (jointly distributed) random variables, which will be important in the coming posts, and really, the most general way we can understand such variables one dimension at a time. We again begin with an approach that seems natural, but doesn't work. For this we just consider two uniform $[0, 1]$ variables and reject if the sum of their squares $> 1$. This should uniformly fill the upper right quarter disk (it's for convenience so that I don't have to deal with a bunch of plus/minus signs with square roots and state all sorts of caveats with inverse trig functions and calculating the integrals. You already know how much I don't like inverse trig functions).

The reasoning for this sampling method is as follows. If we are already given a uniform sample $y$, then the draw that produced it needs to make it fall within the disk, but with no other restrictions. Namely $x$ would be uniformly be within $[0, \sqrt{1-y^2}]$. Thus a seemingly good strategy would be to get the sample $v$, and synthesize $x$ as the appropriate uniform random variable. The density function for $x$ in this case can be computed by some manipulations of indicator functions:

$$\int_0^1 \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) \; dy$$

The reason for this function is that we are simply using integration to take some weighted averages of sample values continuously indexed by $y$. Namely, given a fixed $y$, we have that $x$ is a some possible value for the random variable, with pdf being the conditional distribution of $x$ given $v$. Then the probability density evaluated at $v$ is the weighting in this weighted average. This is called a mixture distribution. By manipulating the inequality $0 \leq x \leq \sqrt{1-y^2}$ to get it in a form that gives its explicit dependence on the variable we actually want to integrate, $y$, we see it is equivalent to $0 \leq y \leq \sqrt{1-x^2}$, not surprising due to the symmetry of a circle:

$$\chi_{[0, \sqrt{1-y^2}]}(x) = \chi_{[0, \sqrt{1-x^2}]}(y).$$

We need this, because a indicator function on $y$ instead of $x$ enables us to just change the limits of integration. Therefore 

$$p_X(x) = \int_0^{\sqrt{1-x^2}} \frac{1}{\sqrt{1-y^2}} \; dy$$

the integrand of which we recognize as the derivative of $\sin^{-1}(y)$. Thus the density is

$$p_X(x) = \sin^{-1}(\sqrt{1-x^2}) - x \sin^{-1}(0) = \sin^{-1}(\sqrt{1-x^2})$$

on $[0, 1]$. What does this actually look like? (To recap: we select the $y$-coordinate with the usual uniform randomness, and then we select the $x$-coordinate uniformly among all the remaining possible values that are possible once $y$ is fixed):

It totally gets scrunched up top there in the upper left. This is because the uniformity in 2 dimensions means that the frequency of sampling of the horizontal coordinate decreases as the vertical coordinate increases, but this is not accounted for if you just sample the vertical coordinate independently first. This is somewhat subtle and not intuitively obvious; it comes from conflicting definitions of uniformity and equality, and that:

One has to realize uniformity does not persist under transformations and projections.

In order to "fix" the problem of the nonuniformity happening here, we consider the uniform distribution in the disk as having density form $\frac{4}{\pi} dx \; dy$. If the conditioning on a specific value of $y$ gives a uniform distribution in $[0, \sqrt{1-y^2}]$, then we have to split it up between the numerator and denominator: essentially all we're doing is "undoing" the extra scrunching. In some sense, it compromises by thinning out what happens vertically, so that things sample less frequently in the small scrunched up horizontal intervals near the upper left corner.

$$\begin{multline}\frac{4}{\pi} dx\; dy = \frac{4}{\pi} \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) dx\; \sqrt{1-y^2}\; dy \\ = \left( \frac{1}{\sqrt{1-y^2}} \chi_{[0, \sqrt{1-y^2}]}(x) dx\right)\left(\frac{4}{\pi} \sqrt{1-y^2}\;dy\right).\end{multline}$$

The first factor is the conditional distribution $p_{X|Y}(x|y)\; dx$, not $p_X(x)\; dx$. Conditional distributions are how you can factor random variable generation into a sequence of generations each depending on the previous result. To evaluate $p_X(x)\; dx$, one integrates the expression over all $y$. This gives $p_X(x)\; dx = \frac{4}{\pi}\sqrt{1-x^2}\; dx$ which is exactly the same a the distribution of $y$! This should not actually be surprising, given that the circle is symmetric between $x$ and $y$ But it should be noted that using this distribution for $x$ after using it for $y$ will not give correct results, simply because once we've simulated $y$, and we act based on that value, the distribution is by definition conditional. How do we simulate $y$ with the distribution $\frac{4}{\pi}\sqrt{1-y^2}$? The (not so fun) cumulative distribution function corresponding to this is:

$$\frac{2}{\pi} \left(y \sqrt{1-y^2} +\sin^{-1}(y)\right)$$ for $y \in [0, 1]$.

We would need to apply the inverse of this function on a uniform random variable to get what we want. This will take a bit of finessing with Newton's method if you really want to solve it. (At the end of this series, I'll have a Jupyter notebook that will allow you to follow along). We've already described how to get the $x$ coordinate from this, but to recap: we just sample uniformly from $0$ to $\sqrt{1-y^2}$ which is simply generating the uniform from $0$ to $1$ and then rescaling.

But the real takeaway here is basically a granular understanding of a formula of conditional probability that is often committed to memory long before building intuition:

$$\rho_{X,Y}(x,y) = \rho_X(x) \rho_{Y|X}(y|x).$$

The first factor is the marginal pdf, namely, what happens when we forget about $Y$ entirely, and the second is the conditional pdf. This formula always holds, and is a replacement of the more classic formula $\rho_{X,Y}(x,y) = \rho_X(x) \rho_Y(y)$ when $X$ and $Y$ are independent (usually learned much earlier). Anyway, in this formula lies the key to bootstrapping ourselves up to higher dimensions! See you in the next post.