r'''Calculate the expectation value of function ``f`` using weighted
    samples (like the output of an importance-sampling run).

    Denoting :math:`x_n` as the sample n and :math:`w_n` as its (normalized)
    weight, the following is returned:

    .. math::

        \sum_{n=1}^{N} w_n f(x_n)
        \mathrm{\ \ where\ \ } \sum_{n=1}^{N}w_n \overset{!}{=} 1

    :param samples:

        Matrix-like numpy array; the samples to be used.

    :param weights:

        Vector-like numpy array; the (unnormalized) importance weights.

    :param f:

        Callable, the function to be evaluated.

def calculate_expectation(samples, weights, f):
    r'''Calculate the expectation value of function ``f`` using weighted
    samples (like the output of an importance-sampling run).

    Denoting :math:`x_n` as the sample n and :math:`w_n` as its (normalized)
    weight, the following is returned:

    .. math::

        \sum_{n=1}^{N} w_n f(x_n)
        \mathrm{\ \ where\ \ } \sum_{n=1}^{N}w_n \overset{!}{=} 1

    :param samples:

        Matrix-like numpy array; the samples to be used.

    :param weights:

        Vector-like numpy array; the (unnormalized) importance weights.

    :param f:

        Callable, the function to be evaluated.

    '''
    assert len(samples) == len(weights), "The number of samples (got %i) must equal the number of weights (got %i)." % (len(samples),len(weights))
    normalization = 0.
    out           = 0.
    for weight, sample in zip(weights, samples):
        normalization += weight
        out += weight * f(sample)
    return out/normalization