def moment(f, pmf, center=0, n=1):
    """
    Return the nth moment of `f` about `center`, distributed by `pmf`.

    Explicitly:   \sum_i (f(i) - center)**n p(i)

    Note, `pmf` is the joint distribution. So n=1 can be used even when
    calculating covariances such as <xx> and <xy>. The first would actually
    be a 2nd moment, while the second would be a mixed 1st moment.

    Parameters
    ----------
    f : array-like
        The numerical values assigned to each outcome of `p`.
    pmf : array-like
        The pmf for a distribution, linear-distributed values.
    center : float
        Calculate a centered moment.
    n : int
        The moment to calculate.

    """
    return ((f - center)**n * pmf).sum()