def update_moments(n, x, m1, m2):
    """Add the n'th sample, x, to the 1st and 2nd moments, m1 and m2.
    E.g. on the first sample, m1 = x1, m2 = 0
    on the second sample, use m1,m2 = update_moments(2, x2, m1, m2)
    m1 is the mean, m2 is the second moment
    sample variance = m2/(n-1), population variance = m2/n
    Parameters:
        n : scalar, sample number (for the first sample, n = 1)
        x : scalar or array_like, sample
        m1 : as x, the previous mean
        m2 : as x, the previous second moment
    Returns:
        m1, m2 : updated mean and second moment

    See:
        Welford, B.P. "Note on a Method for Calculating Corrected Sums of Squares and Products"
        Technometrics Vol. 4, No. 3 (Aug 1962), pp. 419-420
        https://dx.doi.org/10.2307/1266577
    """
    #use Welford's algorithm to maintain stability
    delta = x - m1  # difference from previous mean
    m1 += delta/n
    m2 += delta*(x-m1)
    return m1, m2