Can Two Independent Random Variables Sum to Uniform?

As you well know, I'm a big fan of probability and thinking of it in geometric terms. This is because the sample spaces in probability are almost always some sort of high dimensional state space, that describe the possible outcomes (or configurations) of a system. Probability often boils down to measuring what fraction of the overall total space of outcomes, that your particular experiment resulted in. In all of I've read about and studied over the years, understanding the intuition behind probability has been the theme that recurs the most. Today I want to start on an exploration of the notion of uniformity. Now we actually have visited this theme before as well, and the next couple of posts will actually expand upon those, exploring interesting facts about what happens in higher dimensions. There will be plenty of geometry to visualize!

But today I want to start very simple. What can be simpler than a uniform random variable between 0 and 1? It's basically everyone's random number generator: many distributions are simply calculated by combining multiple uniform random variables in different ways, for example, by arithmetic operations (adding, multiplying, raising to powers), applying functions (exponentials, sin, cos, etc), applying an inverse cumulative distribution function, etc. And we will attempt to do something that may look like it breaks all sorts of rules! Break rules? Consummate rule-follower mathematicians BREAKING RULES? What has the world come to? Now I don't encourage literally breaking rules and deliberately getting into trouble. At least, I'm not supposed to. All I am saying to do is to play around a bit and develop intuition, and try to apply things in ways you're not entirely sure of. It's the sense of adventure and some risk-taking, and what the nature of math exploration actually is. At the end of the day, we will need to have firmer justifications for things and we do need to verify. But, thinking of math as living, breathing, experimental science is not what students are used to thinking. Sometimes the intuitions behind things are lost, as well as a sense of adventure. Let me tell you, this was definitely quite the adventure... You will see me get very excited, and I hope it'll rub off on you, too.

Triangular Distribution

Now, adding independent random variables is a very well-known operation. Given $X$ and $Y$ independent uniform random variables (say on $[0, 1]$), if $Z = X + Y$, this is another random variable. What may be surprising at first is that such a sum is actually NOT uniformly distributed. It is easy enough to see that its maximum range is from 0 to 2. Next, we do something like, ask ourselves, when can the two variables sum to something in the range of 1.5 to 2? Once one of the variables is pinned down, this puts a constraint on the possibilities of the other variable. If we get, say, 0.9 for the first, then the second variable can only be in the range of 0.6 and onward, and actually, you would not actually be able to achieve the value 2, since it requires a value of 1.1, which is not an acceptable value. 

In fact, we get something called the triangular distribution. Numbers close to 1 have the highest probability (density) for a sum, while those closer to the extreme ends have lower probabilities, because they rely on each uniform variable to have values clustered very close to 0 or very close to 1. 

So having been sort of surprised about summing two uniform random variables to be non-uniform, one natural question may be to ask, is it actually possible to sum two independent random variables to be uniform? At least, that's the question that occurred to me. If you like, think about it for a bit and pause... Don't spend more than 10 minutes, though.

....

....

Ok. So it turns out that it actually is impossible. There is a very nice, visual proof at Stats StackExchange, which I actually like a lot. It would have made a nice Nested Tori post by itself. Perhaps I'll try to take it apart and explain it someday, but I also think they give a pretty good exposition as it is.

Is That All There Is to it? Should we just give up?

I mean, yes. The theorem says so. At least, we should give up on actually finding those two random variables... And that is well and good. However, should it be the end of the story? What happens if we try anyway? There's many interesting further explorations one can do here to shore up the probabilistic intuition of what is going on. First, knowing that there ARE random variables that can be subdivided in the manner we just described. The exponential (actually, the negative log of a uniform variable), in particular, is the sum of two identically distributed chi-square random variables. The exponential distribution is also implicated in another somewhat related and very interesting theorem about the product of identically distributed uniform random variables, both raised to a power of another such, discussed here by our Summer of Math host.

A second interesting direction is to realize that there are actually ways to sum random variables to uniform. They just cannot be independent. This has come up in some attempts I made to gain intuition on Bose-Einstein statistics, a quantum phenomenon in which indistinguishable particles or states actually behave as if all that binomial stuff collapses, as if, flipping two coins produced three outcomes of equal probability, {heads, heads}, {heads, tails}, and {tails, tails} and somehow that mixed outcome is weighted less. I could do a whole different post about this, as it is also fun, but fundamentally, it boils down to asking this exact same question about summing to uniform. It's just that, in order to solve it, I had to give up independence, in which I would flip different (unfair) coins depending on the outcome. 

But the third adventure is a real fun one, and that's what the rest of this post is about. Essentially, it had partially been started before learning of that Stats StackExchange proof and I worked it out in different phases. Had I known the StackExchange proof, I might well have never gotten started.

Probability distributions are characterized by their density functions, and we actually very systematically understand what effect does taking a sum of two independent random variables has on their density functions. It's a very well-known operation called convolution.

Convolutions

Given the density functions $\rho_X$ and $\rho_Y$, the density function $\rho_{X+Y}$ is given by the convolution integral

$$\rho_{X+Y}(z) = (\rho_X * \rho_Y)(z) = \int_{-\infty}^\infty \rho_X (x)\rho_Y(z-x) \; dx.$$

The intuition here is using one function to smear out another function by weighted averages. It's somewhat related to the convolutions encountered in machine learning and AI, especially if multidimensional, but convolution is very interesting and important operation in math, useful in all sorts of real analysis contexts. It's also an engineering favorite. This has various nice properties reminiscent of multiplication, like being associative, commutative, and distributive over addition of functions. As for an "identity" ... well... we'll get to that. The reason why convolution is very popular is because of its behavior under the Laplace transform (and Fourier transform as well, but the Laplace transform is what turns out to be the key to really solving this problem). This theorem is the

Convolution Theorem.

Given $f, g$ "nice" functions, the Laplace transform of the convolution, $\mathscr{L}[f*g](s)$, is equal to the product of the Laplace transforms $\mathscr{L}[f](s) \mathscr{L}[g](s)$.

In other words, on the transform side, it goes from being "like" multiplication to just, well, being multiplication. It is for this reason that Laplace transforms are useful in differential equations, and is often quoted as "transforming calculus problems into algebra problems". It should also be noted that convolutions are well-defined for functions that may not necessarily be probability density functions, and the theorem still holds, when the Laplace transforms of such functions exist. Or even, when they don't ... that's a different part of the adventure....

So what is happening in terms of probability density functions (pdfs) is that if one takes Laplace transforms of them, then the Laplace transform (LT) of the pdf of the random variable representing the sum, is the product of the LTs of both pdfs. The LT of uniform with pdf $\chi_{[0,1]}$, the indicator (or characteristic) function that is $1$ on $[0, 1]$ and $0$ elsewhere,

$$\int_{0}^\infty e^{-st} \chi_{[0, 1]} \; dt = \int_{0}^1 e^{-st} \; dt = \frac{1-e^{-s}}{s}.$$

So, now, if we want two independent random variables $\rho_X$ and $\rho_Y$ to represent something that sums to uniform, this means their convolution is the uniform pdf.

But that means ... their Laplace transforms multiply to $ \frac{1-e^{-s}}{s}$! (that's me exclaiming, not taking the factorial). But if they are identically distributed... This means the our purported density function for this is $$\mathscr{L}[\rho_X] = \mathscr{L}[\rho_Y] = \sqrt{\frac{1-e^{-s}}{s}}.$$

Great!!! All we have to do, then, is find the inverse Laplace Transform! No problem at all!!

Except... well... ::opens up a reference book:: looks at the Bromwich Integral... ok where am I supposed to put the contours? What are the singularities? What about that branch cut? Didn't we spend a lot of time last year and years before that hating on branch cuts? ARRRRGH WHAT?? :: slinks away in defeat... ::

Highly Suspect Distributional Calculus

So back to the drawing board. Let's pick apart what that Laplace transform of a uniform pdf really means. It is $\chi_{[0,1]}$, with transform $\frac{1 - e^{-s}}{s}$. As Laplace transform fans know, multiplication by $s$ is differentiation, and division by $s$ is integration starting at some point. Now let's forget that denominator for a moment, and look at pieces separately.

Laplace (and Fourier) transform fans are may be familiar with some (Schwartzian) distributions, or generalized functions (we'll use the latter term since the word "distribution" is overused) namely the Dirac delta (impulse "function") and its derivatives. I actually first learned about that function in an engineering class, not in a math class. And the dang thing is sure useful. The important things to know about this "function" (As a mathematician, I'm obligated to put function in quotes):

1. $\int_{-\infty}^{\infty} \delta(t) \; dt = 1$.
2. $\delta * f = f*\delta = f$ (it is the identity for convolution).
3. The Laplace transform of $\delta$ is $1$ (actually this follows from the above by the Convolution theorem).
4. $\delta$ is the (suitably generalized) derivative of the unit step function (Heaviside function) $\chi_{[0, \infty)}$, the function that is $1$ on the positive real line, and $0$ elsewhere.

It can be visualized as having a graph with infinitely skinny, infinitely high peak at $0$, a value of $0$ elsewhere, yet has a total area under it equal to $1$ (no actual function on the real line can do this: being nonzero at only one point always means it integrates to $0$). When the derivative of something is a $\delta$, that represents a jump discontinuity of size $1$. Yes, I know some of those listed properties do need some more justification, since $\delta$ is not in fact actually a function, but rest assured, there are firm foundations for this.

So if we just look at $1 - e^{-s}$, this corresponds to the Dirac delta "function", minus the Dirac delta "function" shifted to the right by $1$. This is also intuitive: the uniform pdf $\chi_{[0,1]}$  has two jumps, the value going from $0$ to $1$ at the point $0$, and from value $1$ back to $0$, at the point $1$. In general, $e^{-as}\mathscr{L}[f](s)$ the transform of the function $f(x-a)$, another useful property (the translation property)... we'll be using this a lot, soon.

Now if we take the square root:

$$\sqrt{1-e^{-s}}$$

then corresponds to the convolutional square root of that difference of deltas. What can we do with that? There is another thing we know: we know how to expand $\sqrt{1 - e^{-s}}$ in a power series (binomial series) (but crucially, and this is what took some time for me to see, in powers of $e^{-s}$, not $s$). For $\operatorname{Re}(s) > 0$, the series converges uniformly, so there are a number of good analytic operations we can do. Namely, the series is

$$\sqrt{1 - e^{-s}} = (1-e^{-s})^{1/2} = \sum_{n=0}^\infty {1/2 \choose n}(-1)^n e^{-ns}$$

The suspicious looking binomial coefficent ${ 1/2 \choose n}$ is defined using the Gamma function, the usual function that generalizes factorial to non-integers. The notation really is the generalization of binomial coefficient to the non-integer case and is written that way precisely because it generalizes a binomial expresion like $(1+y)^n$ to a non-integer power. It's just that the coefficients, once nonzero, never become zero again, so expanding it out makes it an infinite series. We find that we actually can massage the coefficient a bit to express it in terms of the more normal (integer) binomial coefficients:

$${1/2 \choose n} = \frac{1}{2^{2n-1} n} (-1)^{n+1} { 2n-2 \choose n-1}$$ (and conventionally taking it to be $1$ at $0$). Then $e^{-ns}$ is the transform of $\delta$ translated by $n$ (it's using the property that its LT is 1, and the translation property). Taking this inverse transform we have:

$$\mathscr{L}^{-1}\left[\sqrt{1-e^{-s}}\right] = \sum_{n=0}^\infty { 1/2 \choose n}(-1)^n\delta_n = \delta_0 - \sum_{n=1}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \delta_n =: \rho(t).$$

The latter sum rather nicely gets rid of the signs, so that only the first term has a plus sign and the rest of the terms are subtracted from it.

A more "illegal" way of seeing this is to consider $\delta_0$ the identity operator as a convolution, and $\delta_0 - \delta_1$ as being a "perturbation" of the identity, and we formally take the square root of that difference using Binomial series:

$$\sqrt{\delta_0 - \delta_1} = \delta_0 - \frac{1}{2} \delta_1 - \frac{1}{8} \delta_1^2 - \dots = \sum_{n=0}^\infty { 1/2 \choose n}(-1)^n\delta_1^n$$

and using the fact that $\delta_1^n = \delta_n$ (because convolution by $\delta_1$ multiple times will translate by 1 over and over again). This is the best of Euler and he would do it with abandon, without regard to convergence. The advantage of the Laplace transform method is that we have all the series converging and we know that stuff like term-by-term integration is guaranteed to work in such circumstances.

In either case, this is indeed a generalized function in the same sense as $\delta$ is, and in fact, a tempered one: given a function $\phi$ the very nice Schwartz class of smooth functions vanishing at infinity after multiplication by polynomial of any degree, when $\rho$ operates on it, the resultant infinite series of translations of $\phi$ must converge, since each discrete value will be dominated by any negative $n$th power. The situation is even nicer with smooth functions of compact support (vanishing outside a finite interval): $\rho$ will always be a finite sum at any evaluation point. Suffice to say, the convergence is good enough.

Next We Will Need: The Square Root of .... Integral?

Ah. But we're missing something. We said to ignore that denominator $s$, so we're still not done yet. It's time to bring it back. Now, just like multiplying by $s$ corresponds to differentiation, dividing by $s$ corresponds to integration (from 0). Since the convergence is nice, dividing by $s$ we can just integrate that crazy sum above term-by-term, and be done..........

Oh wait. We divide by the square root of $s$, not $s$. Oops. What does dividing by the square root of $s$ mean? Dividing by $s^2$ is integrating twice. So dividing by the square root of $s$ is, of course, just integrating half a time. WHAT??? Yes the square root of integral!! (a double exclamation. Not the double factorial of integral, although I wonder if it's possible to..........)

One additional way to motivate things is that differentiation plays very well with convolution. Namely $\partial^{i+j} ( f * g) = (\partial^{i+j} f) * g = f*(\partial^{i+j} g) = \partial^i f * \partial^j g$ in any combination that makes sense. You can split the differentiations up to happen on each component that can take that many derivatives, you can put them all on one factor, or the other factor, you can re-split it up in different ways, etc. Similarly, if we think of integration from $0$ to $t$ as a negative order differential, it all continues to work out, too.

So here if we define $\rho$ to be the binomial-series sum of deltas above, we have $\rho * \rho = \delta_0 - \delta_1$ by construction. If we have to do it on a test function, fine. The second convolution will stagger the additional terms, so they all cancel exactly outside of $[0, 1]$. Then the function we seek, $\tilde \rho$ satisfies $\tilde \rho * \tilde \rho = \chi_{[0,1]} = \int_0^x \delta_0(t) - \delta_1(t) \; dt$.

Thus $\tilde \rho * \tilde \rho = \int_0^x \rho * \rho \; dt = \partial^{-1}(\rho * \rho)$, that is, it is the integral of the whole thing. For our next illegal trick, we're going to split that integration in exactly half, and put it on each individual $\rho$:

$$\tilde \rho * \tilde \rho = \partial^{-1}(\rho * \rho) = \partial^{-1/2} \rho * \partial^{-1/2} \rho.$$

Total nonsense?? Maybe. I have not actually justified for generalized functions in general, but we'll see in a short bit how Laplace transforms take care of all of it, at least with the exponent of ½. For functions vanishing at $0$, it sure is convenient (vanishing at $0$ makes it work well with Laplace transforms, since we often need to subtract derivatives of initial conditions). It turns out the iterated integral operation has, what is called, a kernel representation: you can represent iterated integrals as one single integral. The "kernel" in kernel representation is analogous to continuous version of a matrix, and the integral against a function is like a continuous sum version of matrix multiplication.

The kernel representation of iterated integration is:

$$\int_0^t \int_0^{t_{n-1}} \dots \int_0^{t_0} f(t_0) \; dt_0 \; dt_1 \; \dots \; dt_{n-1} = \frac{1}{(n-1)!}\int_{0}^t f(\tau)(t-\tau)^{n-1}\; d\tau.$$

The proof is actually very pretty straightforward: interchanging the order of integration and proof by induction. When the innermost variables are switched, the inside is constant with respect to the switched variable, so integrates to a difference. Stacking them up will give you the powers, so now it's only one integral.

But the amazing thing is... that expression is defined for all real numbers $\nu$ (we can rely on our trusty old friend, the Gamma function, for the factorial). We DEFINE the fractional integral

$$I^{\nu}(f)(t) = \partial^{-\nu}(f)(t) = \frac{1}{\Gamma(\nu)}\int_0^t (t - \tau)^{\nu - 1}f(\tau) \; d\tau.$$

The important thing to know here is $I^{\nu + \mu} = I^\nu I^\mu$. This can be shown by composing two integrals and using beta functions (of Bayesian statistics fame) to rewrite the product $\Gamma(\alpha) \Gamma(\beta)$. It will work for non-integer $\nu$ and $\mu$. It's a bit of a mess, so we won't reproduce all of it here. Anyway, this makes integration operations of all orders a form of convolution—the telltale sign being, the $\tau$ and the $t-\tau$, though there is some funny business with the limits of integration. Also, there's trouble if $\nu$ is a nonpositive integer (because the gamma function is infinite in such places), with the function $\frac{1}{\Gamma(\nu)}t^{\nu-1} \chi_{[0, \infty)}(t)$. Thus if we want to write it for $\nu = -\frac{1}{2}$, we have (recalling $\Gamma(1/2) = \sqrt{\pi}$, another amazing fact, and another place where $\pi$ has a habit of just inviting itself to random parties):

$$I^{1/2}(f)(t) = \frac{1}{\sqrt{\pi t}} \chi_{[0, \infty)}(t) * f(t).$$

We can also justify this a bit more by showing that the Laplace transform of $\frac{1}{\sqrt{\pi t}}\chi_{[0, \infty)}(t)$ is in fact $\frac{1}{\sqrt{s}}$ so that it just becomes $\frac{1}{\sqrt{s}}$ times the original series with $\sqrt{1-e^{-s}}$. Then this means the original transform is just a series of translated $\frac{1}{\sqrt{\pi t}}\chi_{[0, \infty)}(t)$'s. But STILL, ISN'T IT TOTALLY AWESOME THAT WE TOOK THE SQUARE ROOT OF INTEGRAL???

AND FINALLY...

So what's left to do? We need to take the ½th integral of that sum of delta's above. So this is just convolving the function $$\frac{1}{\sqrt{\pi t}} \chi_{[0, \infty)}(t)$$ with the majestic

$$\delta_0 - \sum_{n=1}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \delta_n$$

Convolutions with deltas are supposed to be trivial, because they're the identity op, remember? We have some shifted deltas, however. Convolving with it gives a shifted function, namely, we shift that $\frac{1}{\sqrt{\pi t}} \chi_{[0,\infty)}(t)$ over by $n$ when we apply $\delta_n$. The result is therefore

$$\tilde \rho(t) = \frac{1}{\sqrt{\pi t} } \chi_{[0, \infty)}(t)* \rho(t) = \frac{1}{\sqrt{\pi t} }\chi_{[0, \infty)}(t)- \sum_{n=1}^\infty\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \frac{1}{\sqrt{\pi (t-n)} }\chi_{[0, \infty)}(t-n)$$

Those indicator functions are bloody annoying looking. I kept them in specifically because it allows you to keep track of the bounds of integration better. But there's a nice reformulation, since here, the action of these indicator functions is merely to cut off the sums depending on whether we have encountered the relevant delta term. In fact, its effect is to finitize the sum (which turns out to allow us to dodge some convergence questions): at any given $t$, only the terms up to $\lfloor t \rfloor$ have been added! We thus rewrite:

$$\tilde \rho(t) = \frac{1}{\sqrt{\pi t} } - \sum_{n=1}^{\lfloor t \rfloor}\frac{1}{2^{2n-1} n}{ 2n-2 \choose n-1} \frac{1}{\sqrt{\pi (t-n)} }$$

for $t > 0$ and $0$ otherwise. This includes a value of $\pm \infty$ at every integer. To avoid that, we can use $\lceil t \rceil -1$ instead of $\lfloor t \rfloor$, but infinite values on a discrete set of measure zero (set with no length) don't affect the result, especially when it comes to convolution, so we'll not bother.

Of course, we should actually probably take a look at what we just constructed. This is the graph of $\tilde{\rho}$:

You can very plainly see it has negative values, which prevents it from being a pdf. We've shown that a convolutional square root does in fact exist, but it is not a pdf. Because the uniform pdf has a downward jump, no fix can be made to eliminate the negative values (since a downward jump requires convolution by something with vertical asymptotes, whose values approach $-\infty$). And that, my friends, is my very satisfying negative answer to the this post's title! The the number of different tools and areas of math that this touches upon, not to mention, having just a random factor of $\pi$ showing up where there are no circles, is a real delight!