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

Question

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?

Answer 1

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
        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()