How does random() * random() distribution look like?

Question

Given a function random() which returns floating-point value uniformly distributed between 0 and 1.

What's the type of distribution of function random() * random()?

Answer 1

# test.py
import numpy as np
import matplotlib.pyplot as plt

N = 10**6
plt.hist(np.random.uniform(size=N) * np.random.uniform(size=N), bins=50, normed=True)
plt.show()

Running python test.py produces:

Answer 2

The type is a Product Distribution, it is not uniform anymore.

Answer 3

The transformation is y = x*x. x has the probability distribution function fx(x) = 1 on the range 0 <=x <= 1. The cumulative distribution function of x is then Fx(x) = x on the same range.

The CDF of y is Fy(y <= Y) = Fx(sqrt(Y)) = sqrt(Y), 0 <= Y <= 1.

Now differentiate to get fy(y) = 1/(2*sqrt(y)) on the same range.

EDIT:

The above solution assumes that "random() * random()" uses the same value dependently for each draw. If instead you want the multiplied values to be independent of one another, the math is more involved but still tractable.

Now let

y1 = x1*x2 where fx1(x1) = 1 on 0 <= x1 <= 1 and similarly for x2.

Assuming independence between x1 and x2, the joint PDF is

fx1x2(x1,x2) = fx1(x1)*fx2(x2).

Introduce an additional variable y2 to deal with the transformation of 2-variable joint PDFs. For nice calculations, let y2 = x2.

So our system is

g1(x1,x2) = x1*x2

g2(x1,x2) = x2

As in the simpler case, we need to invert the function, now by solving for both y1 and y2:

h2(y1,y2) = x2 (= y2)

h1(y1,y2) = y1/x2 = y1/y2

We will need the Jacobian

J = (pg1/px1)(pg2/px2) - (pg1/px2)(pg2/x1)

where "p" is partial derivative.

So in our case

J = (x2)(1) - (x1)(0) = x2.

The transformation formula (from any intro calculus-based probability text) is

fy1y2(y1,y2) = fx1x2(x1,x2)/J

which in our case simplifies to

1/y2 on the range 0 <= y1/y2 <= 1 and 0 <= y2 <= 1.

Finally, to get fy1(y1) we integrate the joint distribution over the unwanted variable y2, taking care to stay in the correct range y1/y2 <= 1 or y1 <= y2 since y2 >= 0.

fy1(y1) = Integral from y1 to 1 of (1/y2)dy2 = -ln(y1) on the range 0 <= y1 <= 1.

Note that in both cases, the distribution of the product is weighted in favor of smaller values, since a fraction (0 <= x <= 1) times a fraction is a smaller fraction.