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
        cdf_med = scipy.integrate.quad(fun, minus_inf, med)[0]
        if cdf_med < q:
            lower = med
            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 """
    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:
    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')

enter image description here

  • $\begingroup$ Interesting. So, the "right" approach to compute statistics for a PDF is to draw samples from it, rather that directly evaluating it? $\endgroup$ – Konstantin Jan 15 '18 at 10:43
  • $\begingroup$ Direct evaluation is usually faster when you can do, but here it is difficult. $\endgroup$ – Björn Jan 15 '18 at 11:25
  • $\begingroup$ @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. $\endgroup$ – David Dale Jan 15 '18 at 11:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.