Skip to content

Instantly share code, notes, and snippets.

@lzkelley
Created May 11, 2019 01:21
Show Gist options
  • Save lzkelley/03443273bf8a867716a9247ce01d49f5 to your computer and use it in GitHub Desktop.
Save lzkelley/03443273bf8a867716a9247ce01d49f5 to your computer and use it in GitHub Desktop.
import numpy as np
import scipy as sp
import scipy.stats
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=[8, 5])
ax.set_xscale('log')
# Construct a log-normal distribution
dist = sp.stats.lognorm(1)
# Samples from the distribution
x = dist.rvs(1000)
# fine-spacing evaluation points
xs = np.logspace(-2, 2, 2000)
# coarse-spacing bins
bins = np.logspace(-2, 2, 40)
# Plot "true" PDF and binned samples in (BLACK)
ax.hist(x, bins=bins, color='k', rwidth=0.9, alpha=0.5, density=True)
ax.plot(xs, dist.pdf(xs), 'k-', lw=2)
# Construct normal, linear, gaussian KDE
kde = sp.stats.gaussian_kde(x)
# Plot KDE PDF at evaluation points (RED)
yy = kde(xs)
ax.plot(xs, yy/np.max(yy), 'r-', lw=2.0)
# Draw from KDE and plot histogram
bb = kde.resample(1000)[0]
ax.hist(bb, bins=bins, color='r', rwidth=0.5, alpha=0.5, density=True)
# Construct a KDE in log-space
sp_kde_log = sp.stats.gaussian_kde(np.log10(x))
# Get PDF, multiply by inverse of Jacobian (1/x)
yy = sp_kde_log(xs)*xs
# Plot log-space KDE PDF (BLUE)
ax.plot(xs, yy/np.max(yy), 'b-', lw=2)
# Sample (by inverting cumulative distribution)
zz = np.cumsum(yy)/np.sum(yy)
aa = np.interp(np.random.uniform(0.0, 1.0, 1000), zz, xs)
# Plot histogram of log-space KDE drawn samples (BLUE)
ax.hist(aa, bins=bins, color='b', rwidth=0.5, alpha=0.5, density=True)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment