单遍通用统计的整数数值稳定性：多统计量单遍算法疑问

Great question—numerical stability is one of those easy-to-overlook details that can totally derail your stats calculations, especially with large datasets or values clustered tightly around a large mean. Let’s break this down step by step:

First: Why the Berkeley Method Fails #

The formula you mentioned for standard deviation (std = sqrt((sum(x²)-N*mean(X)^2)/(N-1))) works mathematically, but it’s prone to catastrophic cancellation—the main culprit behind numerical instability here.

Here’s what happens: When your dataset’s mean is very large, sum(x²) and N*mean(X)^2 become two extremely large, almost identical numbers. Subtracting them erases most of the significant digits in the result. For example:

• Suppose every data point is 1e6 + 0.1 (so mean is 1e6 + 0.1)
• sum(x²) for N=1000 is 1000*(1e6+0.1)^2 = 1e15 + 2e7 + 10
• N*mean(X)^2 is exactly the same value (since all points are identical)
• But due to floating-point precision limits, these two values might be stored as slightly rounded versions. Subtracting them could give you 0 instead of the correct 10, leading to a standard deviation of 0—completely wrong.

This issue gets worse as the mean grows relative to the spread of your data.

Stable One-Pass Alternatives #

The gold standard for stable one-pass calculation of these stats is extensions of Welford’s Online Algorithm, which avoids catastrophic cancellation by updating stats incrementally using small, per-point differences instead of large cumulative sums.

1. Mean & Standard Deviation #

Welford’s algorithm tracks the current mean and a cumulative term M2 (which represents the sum of squared deviations from the mean, updated incrementally):

def welford_update(n, mean, M2, x):
    n += 1
    delta = x - mean
    mean += delta / n
    delta2 = x - mean
    M2 += delta * delta2
    return n, mean, M2

# After processing all points:
std_dev = sqrt(M2 / (n - 1)) if n > 1 else 0.0

• Why this works: Each update uses delta (the difference between the new point and the current mean), which is a small number relative to x or mean. No large numbers are subtracted, so precision is preserved.

2. Skewness & Kurtosis #

You can extend Welford’s algorithm to track third (M3) and fourth (M4) moments for skewness and kurtosis. Here’s the incremental update logic:

def update_skew_kurt(n, mean, M2, M3, M4, x):
    n += 1
    delta = x - mean
    delta_n = delta / n
    delta_n2 = delta_n * delta_n
    term1 = delta * delta_n * (n - 1)
    
    # Update mean first
    mean += delta_n
    
    # Update higher moments
    M4 += term1 * delta_n2 * (n*n - 3*n + 3) + 6*delta_n2*M2 - 4*delta_n*M3
    M3 += term1 * delta_n * (n - 2) - 3*delta_n*M2
    M2 += term1
    
    return n, mean, M2, M3, M4

# Final calculations:
variance = M2 / n
skewness = (M3 / n) / (variance ** (3/2)) if variance != 0 else 0.0
# Excess kurtosis (subtract 3 to get relative to normal distribution)
kurtosis = (M4 / n) / (variance ** 2) - 3 if variance != 0 else 0.0

• Note: Skewness and kurtosis definitions can vary slightly (e.g., sample vs population), so adjust the divisors (n vs n-1) based on your needs.

3. Covariance #

For two variables x and y, use a similar incremental approach to track their means and a cumulative covariance term Cxy:

def update_covariance(n, mean_x, mean_y, Cxy, x, y):
    n += 1
    delta_x = x - mean_x
    delta_y = y - mean_y
    mean_x += delta_x / n
    mean_y += delta_y / n
    Cxy += delta_x * (y - mean_y)
    return n, mean_x, mean_y, Cxy

# Final covariance:
covariance = Cxy / (n - 1) if n > 1 else 0.0

• This avoids the unstable sum(x*y) - N*mean_x*mean_y formula by using incremental differences, just like Welford’s method.

Key Takeaway #

All these stable algorithms work by never subtracting two large numbers—instead, they build up stats using small, per-point deviations from the current mean. This preserves floating-point precision even for large datasets or means, eliminating the numerical instability you’d get with the Berkeley method.

内容的提问来源于stack exchange，提问作者zlon