import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
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(1000000)])
samples.sort()
print 'First 5', samples[:5]
print 'Last 5', samples[-5:]
print 'Sample Median', np.median(samples)
First 5 [-190223.22139462 -181277.24263936 -147879.31383188 -74596.48367245 -31040.63070675] Last 5 [ 18013.45871543 18739.73304085 19236.4912174 20586.83355282 103165.42384118] Sample Median -1.99954147793
fig, ax = plt.subplots(1, 2, figsize=(12,4))
ax[0].hist(samples, bins=100, range=(-6, 6))
ax[0].set_title('Histogram of my random deviates')
scipy_cauchy_pdf = scipy.stats.cauchy(loc=cauchy_mu).pdf
xs = np.linspace(-6, 6, num=1000)
ax[1].plot(xs, scipy_cauchy_pdf(xs))
ax[1].set_title('Scipy PDF')
<matplotlib.text.Text at 0xb7486ec>