I show that the log-normal distribution converges to 1/x as sigma approaches infinity by plotting 1/x along with the log-normal distribution for increasing values of sigma. Since the log-normal distribution is normalized to sum to one but 1/x is not, I plot the log-normal distributions unnormalized so they can be visually compared to 1/x.
import numpy as np
import matplotlib.pyplot as plt
def unnormalized_lognormal_pdf(x, mu, sigma):
"""Unnormalized probability density function for the log-normal distribution.
The 'np.sqrt(2*np.pi) * sigma' has been removed from the denominator, so the
result is no longer a valid probability distribution because it does not sum
to one. The shape is still the same though because the missing normalizer
was just a constant wrt x.
"""
normalizer = 1.0 / x
exp_part = np.exp(-((np.log(x) - mu)**2)/(2 * sigma**2))
return normalizer * exp_part
vectorized_lognormal_pdf = np.vectorize(unnormalized_lognormal_pdf)
def plot_unnormalized_lognormal(x, mu, sigma):
"""Plots the unnormalized log-normal distribution with the specified parameters."""
label = "unnormalized lognormal(mu={0}, sigma={1})".format(mu, sigma)
plt.plot(x, vectorized_lognormal_pdf(x, mu, sigma), label=label)
NUM_X_POINTS = 10000
# Plot 1/x along with unnormalized lognormals with increasing sigma.
x = np.linspace(1e-6, 1, num=NUM_X_POINTS)
plt.loglog()
one_over_x = 1.0 / x
plt.plot(x, one_over_x, label='1/x')
for sigma in (2, 3, 4, 8, 16):
plot_unnormalized_lognormal(x, 0, sigma)
plt.legend(loc='lower right')
plt.xlabel('x')
plt.ylabel('p(x)')
<matplotlib.text.Text at 0x973176c>
Another way to show the result is to calculate the mean absolute difference between 1/x and the unnormalized log-normal distribution at several points between 0 and 1. The difference approaches zero as sigma increases.
sigmas = []
mean_abs_errors = []
for exponent in xrange(0, 20):
sigmas.append(2**exponent)
# Average the difference between 1/x and the lognormal at `NUM_X_POINTS` points.
mean_abs_errors.append(np.mean(np.abs(one_over_x - vectorized_lognormal_pdf(x, 0, 2**exponent))))
for sigma, mean_abs_error in zip(sigmas, mean_abs_errors):
print "sigma {0}, error {1}".format(sigma, mean_abs_error)
plt.loglog()
plt.plot(sigmas, mean_abs_errors)
plt.xlabel("sigma")
plt.ylabel("mean absolute error from 1/x")
-c:10: RuntimeWarning: overflow encountered in long_scalars -c:10: RuntimeWarning: divide by zero encountered in double_scalars
sigma 1, error 108.516971028 sigma 2, error 107.263781685 sigma 4, error 104.57058235 sigma 8, error 79.4518416827 sigma 16, error 31.6909090087 sigma 32, error 9.04754313241 sigma 64, error 2.34061429391 sigma 128, error 0.590218129844 sigma 256, error 0.147873369473 sigma 512, error 0.0369883058605 sigma 1024, error 0.00924832474948 sigma 2048, error 0.00231215921399 sigma 4096, error 0.000578044680308 sigma 8192, error 0.00014451147488 sigma 16384, error 3.6127887772e-05 sigma 32768, error 9.03197313323e-06 sigma 65536, error 2.25799335929e-06 sigma 131072, error 5.64498345743e-07 sigma 262144, error 1.41124589157e-07 sigma 524288, error 3.52811502797e-08
<matplotlib.text.Text at 0x9f19d6c>