import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
np.random.seed(42)
def cauchy_inverse_cdf(p, mu, sigma):
return mu + sigma * np.tan(np.pi * (p - 0.5))
def get_cauchy_deviate(mu, sigma):
"""Get a random cauchy deviate using the transformation method."""
uniform_deviate = np.random.random()
return cauchy_inverse_cdf(uniform_deviate, mu, sigma)
cauchy_mu = -2
cauchy_sigma = 0.5
samples = np.array([get_cauchy_deviate(cauchy_mu, cauchy_sigma) for _ in xrange(10000000)])
samples.sort()
print 'First 5', samples[:5]
print 'Last 5', samples[-5:]
print 'Sample Median', np.median(samples)
First 5 [-472102.55301209 -421651.84579584 -415302.55927885 -411290.46997529 -407346.13437632] Last 5 [ 340251.58424335 551800.27483741 701783.98767175 1178281.1345642 1906305.9508952 ] Sample Median -1.99996602588
fig, ax = plt.subplots(figsize=(8, 8))
ax.hist(samples, bins=100, range=(-6, 2), label='Deviates', normed=True, alpha=0.5)
scipy_cauchy_pdf = scipy.stats.cauchy(loc=cauchy_mu, scale=cauchy_sigma).pdf
xs = np.linspace(-6, 2, num=1000)
ax.plot(xs, scipy_cauchy_pdf(xs), color='red', label='Scipy PDF', linewidth=3)
ax.legend()
<matplotlib.legend.Legend at 0xa9e3dcc>