def calculate_expectation(samples, weights, f):
    """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.0
    out = 0.0
    for weight, sample in zip(weights, samples):
        normalization += weight
        out += weight * f(sample)
    return out / normalization
