import matplotlib.pyplot as plt
import numpy as np
np.random.seed(42)
Estimate and plot the number of U(0,1) numbers you can multiply before underflowing.
num_prods = [np.nonzero(np.cumproduct(np.random.random(1000)) == 0)[0][0] + 1 for _ in xrange(200000)]
mean = np.mean(num_prods)
sd = np.std(num_prods)
print "{0} +/- {1} products before underflowing".format(mean, sd)
rounded_mean = np.round(mean)
plt.hist(num_prods, bins=np.arange(rounded_mean - 3*sd, rounded_mean + 3*sd, 4))
plt.title("Histogram of the Number of Products Before Underflow")
746.410685 +/- 27.3135432127 products before underflowing
<matplotlib.text.Text at 0xaf64d4c>
I first estimate the mean and variance of log(U). It seems that the mean is -1 and the variance is -1.
log_uniform = np.log(np.random.random(1000000))
mu = np.mean(log_uniform)
var = np.var(log_uniform)
mu, var
(-1.0005811922459742, 1.0003386611570766)
By the central limit theorem, the sum of N random variables $log(U)$ should be distributed as $\mathcal{N}(-N, N)$. Experimentally, this seems accurate.
num_sums = 10000
sum_logs = [np.sum(np.log(np.random.random(num_sums))) for _ in xrange(20000)]
mean = np.mean(sum_logs)
sd = np.std(sum_logs)
print "{0} +/- {1} (variance {2})".format(mean, sd, sd**2)
plt.hist(sum_logs, bins=25)
plt.title('Histogram of sum of {0} log(U)'.format(num_sums))
-9998.68065929 +/- 100.524298534 (variance 10105.1345958)
<matplotlib.text.Text at 0xb0ba34c>
In part (a), I estimated that you can multiply an average of about 746 $U(0,1)$ values before underflowing.
As we know, adding the log of numbers and then exp'ing them is the same as multiplying the numbers. In part (b), I showed that the mean of the sum of $N$ $log(U)$ values is $-N$. Thus we would expect that we could sum about 746 $log(U)$ values before the exp'ing would underflow
for s in xrange(-740, -750, -1):
print s, np.exp(s)
-740 4.19955798965e-322 -741 1.53160350211e-322 -742 5.43472210425e-323 -743 1.97626258336e-323 -744 9.88131291682e-324 -745 4.94065645841e-324 -746 0.0 -747 0.0 -748 0.0 -749 0.0
WHAMMY!