<jupyter_start><jupyter_text>Density estimation demoHere we demonstrate how to use the ``inference.pdf`` module for estimating univariate probability density functions from sample data.<jupyter_code>from numpy import linspace, zeros, exp, log, sqrt, pi
from numpy.random import normal, exponential
from scipy.special import erfc
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Kernel-density estimationGaussian kernel-density estimation is implemented via the `GaussianKDE` class:<jupyter_code># generate some sample data to use as a test-case
N = 150000
sample = zeros(N)
sample[:N//3] = normal(size=N//3)*0.5 + 1.8
sample[N//3:] = normal(size=2*(N//3))*0.5 + 3.5
# GaussianKDE takes an array of sample values as its only argument
from inference.pdf import GaussianKDE
PDF = GaussianKDE(sample)<jupyter_output><empty_output><jupyter_text>Instances of density estimator classes like `GaussianKDE` can be called as functions to return the estimate of the PDF at given spatial points:<jupyter_code>x = linspace(0, 6, 1000) # make an axis on which to evaluate the PDF estimate
p = PDF(x) # call the instance to get the estimate<jupyter_output><empty_output><jupyter_text>We could plot the estimate manually, but for convenience the `plot_summary()` method will generate a plot automatically as well as summary statistics:<jupyter_code>PDF.plot_summary()<jupyter_output><empty_output><jupyter_text>The summary statistics can be accessed via properties or methods:<jupyter_code># the location of the mode is a property
mode = PDF.mode

# The highest-density interval for any fraction of total probability is returned by the interval() method
hdi_95 = PDF.interval(frac = 0.95)

# the mean, variance, skewness and excess kurtosis are returned by the moments() method:
mean, variance, skewness, kurtosis = PDF.moments()<jupyter_output><empty_output><jupyter_text>By default, `GaussianKDE` uses a simple but easy to compute estimate of the bandwidth (the standard deviation of each Gaussian kernel).However, when estimating strongly non-normal distributions, this simple approach will over-estimate required bandwidth.In these cases, the cross-validation bandwidth selector can be used to obtain better results, but with higher computational cost.<jupyter_code># to demonstrate, lets create a new sample:
N = 30000
sample = zeros(N)
sample[:N//3] = normal(size=N//3)
sample[N//3:] = normal(size=2*(N//3)) + 10

# now construct estimators using the simple and cross-validation estimators
pdf_simple = GaussianKDE(sample)
pdf_crossval = GaussianKDE(sample, cross_validation = True)

# now build an axis on which to evaluate the estimates
x = linspace(-4,14,500)

# for comparison also compute the real distribution
exact = (exp(-0.5*x**2)/3 + 2*exp(-0.5*(x-10)**2)/3)/sqrt(2*pi)

# plot everything together
plt.plot(x, pdf_simple(x), label = 'simple')
plt.plot(x, pdf_crossval(x), label = 'cross-validation')
plt.plot(x, exact, label = 'exact')
plt.ylabel('probability density')
plt.xlabel('x')
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Functional density estimation for unimodal PDFsIf we know that the distribution being estimated is a single (but potentially highly skewed) peak, the `UnimodalPdf` class can robustly estimate the PDF even at smaller sample sizes. It works by fitting a heavily modified Student-t distribution to the sample data.<jupyter_code># Create some samples from the exponentially-modified Gaussian distribution
L = 0.3 # decay constant of the exponential distribution
sample = normal(size = 3000) + exponential(scale = 1./L, size = 3000)

# create an instance of the density estimator
from inference.pdf import UnimodalPdf
PDF = UnimodalPdf(sample)

# plot the estimate along with the exact PDF for comparison
x = linspace(-5, 15, 1000)
exact = 0.5*L*exp(0.5*L*(L-2*x))*erfc((L-x)/sqrt(2)) # exact PDF for the exp-gaussian distribution

plt.plot(x, PDF(x), label = 'UnimodalPdf estimate', lw = 3)
plt.plot(x, exact, label = 'exact distribution', ls = 'dashed', lw = 3)
plt.ylabel('probability density')
plt.xlabel('x')
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()<jupyter_output><empty_output>