docs: 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