Compute parameters of a PDF (probability density function) for which no closed form expression is available

I would like to compute parameters such as mean, variance, quantiles, etc. for a PDF which is only given as a piece of code. That is, it can only be evaluated numerically at given points; no closed-form expression.

For example, after using scikit-learn for kernel density estimation, I would like to compute parameters of the resulting PDF.

I feel that popular libraries such as numpy/scipy, scikit-learn, or pandas probably provide such functionality. However, I couldn't find it.

Can anybody please point me to it?

Your problem could be solved either by direct numeric integration or by MCMC.

Numeric integration can be performed most easily by scipy:

import numpy as np
import scipy.stats
import scipy.integrate

def weird_density(x):
""" The function I want to sample """
return scipy.stats.lognorm.pdf(x, 1.0)

def quad_quantile(fun, q, precision=1e-10, minus_inf=-10e10, lower=-10e10, upper=10e10):
""" Use bisection to evaluate a quantile """
while upper - lower > precision:
med = (upper + lower) * 0.5
if cdf_med < q:
lower = med
else:
upper = med
return (upper + lower) * 0.5

print('true mean  :', scipy.stats.lognorm(1).mean())
print('integrated mean:', scipy.integrate.quad(lambda x: weird_density(x) * x, 0, 100)[0])
print('true median  :', scipy.stats.lognorm(1).ppf(0.5))
print('integrated median:', quad_quantile(weird_density, 0.5, 0.001, 0, 0, 100))


The output is

true mean  : 1.6487212707
integrated mean: 1.6484641126903046
true median  : 1.0
integrated median: 0.9998321533203125


However, if the distribution is multidimensional, direct integration may be intractable. In this case, MCMC methods (Markov chain Monte Carlo) can help.

What they do is just make a correct (in some sense) sample from the distribution for which you know the PDF. When you have a sample, you can calculate all your parameters from it as classical sample statistics, just like from any observed data.

MCMC algorithms may be implemented manually, like in the example below. Another option is to use the PyMC3 library or its analogs.

One of the best known MCMC algorithms, Metropolis-Hastings, works as follows:

def sample_next(x):
""" Generate random next point"""
return scipy.stats.norm(loc=x).rvs(1)[0]

def metropolis_hastings(density, generator, n_samples, starting_point=1, random_state=None):
""" Generate sample from the given density function """
np.random.seed(random_state)
result = [starting_point]
for i in range(n_samples):
current_point = result[-1]
next_candidate = generator(current_point)
acceptance_ratio = density(next_candidate) / density(current_point)
if np.random.uniform() <= acceptance_ratio:
result.append(next_candidate)
else:
result.append(current_point)
return np.array(result)

sample = metropolis_hastings(weird_density, sample_next, 1000, random_state=1)

print('sample mean:', np.mean(sample))
print('sample median:', np.median(sample))


This code prints

sample mean: 1.2872857165
sample median: 0.932826883898


You see that the results are not very close, but with longer sampling (1000 points is too few) they will converge. And the shape of distribution is already well matched:

import matplotlib.pyplot as plt
plt.hist(sample, bins=30, normed=True)
idx = np.linspace(0, sample.max())
plt.plot(idx, scipy.stats.lognorm(1.0).pdf(idx))
plt.legend(['true density', 'sample density'])
plt.title('convergence of MCMC sample')
plt.show()


• Interesting. So, the "right" approach to compute statistics for a PDF is to draw samples from it, rather that directly evaluating it? Jan 15, 2018 at 10:43
• Direct evaluation is usually faster when you can do, but here it is difficult. Jan 15, 2018 at 11:25
• @Konstantin, in one-dimensional case you can find quantiles and moments by integration. I have updated the answer. But in multidimensional case, integration is usually intractable, so sampling is the only option. Jan 15, 2018 at 11:53