有没有比sum（）/N更好的算术平均数计算方法？

Question

假设我们有N个数字（整数，浮点数或其他任何类型），并且想要找到它们的算术平均值。最简单的方法是将所有值相加，然后除以值的数量：

def simple_mean(array[N]): # pseudocode
    sum = 0
    for i = 1 to N
       sum += array[i]
    return sum / N

这段代码运行良好，但需要使用大整数。如果我们不想使用大整数，并且可以容忍四舍五入误差，并且N是2的幂次方，我们可以使用“分而治之”的方法：((a+b)/2 + (c+d)/2)/2 = (a+b+c+d)/4，((a+b+c+d)/4 + (e+f+g+h)/4)/2 = (a+b+c+d+e+f+g+h)/8，以此类推。

def bisection_average(array[N]):
   if N == 1: return array[1]
   return (bisection_average(array[:N/2])+bisection_average(array[N/2:]))/2

还有其他的方法吗？

附注：懒人专用平台

- maykeye

---

有趣，但是“可以容忍舍入误差”的那一部分让我感到担忧。我更喜欢没有误差的方法。 - pavium

转念一想，我明天再回来看看，如果我仍然确信它没有严重错误，那么我会撤销删除我的答案... - Matthew Scharley

@pavium：如果你想要一个没有错误的方法，你需要手动计算。 - MusiGenesis

@MusiGenesis -- 嘿，如果我自己手算这个，那么会保证有K个错误，其中K >> 1，而K/N可能逐渐接近某个常数，如1/e之类的东西。;) - Jason S

Answer 1

Knuth列出了以下方法，用于计算给定浮点数的平均值和标准差（原文在1998年版计算机程序设计艺术第2卷的第232页上；我下面的改编避免了对第一次迭代进行特殊处理）：

double M=0, S=0;

for (int i = 0; i < N; ++i)
{
    double Mprev = M;
    M += (x[i] - M)/(i+1);
    S += (x[i] - M)*(x[i] - Mprev);
}

// mean = M
// std dev = sqrt(S/N) or sqrt(S/N+1)
// depending on whether you want population or sample std dev

- Jason S

---

“S += (x[i] - M)(x[i] - Mprev);” 应该改为 “S += (x[i] - Mprev)(x[i] - Mprev);” 吗？ - Gábor Bakos

不。请参阅http://jonisalonen.com/2013/deriving-welfords-method-for-computing-variance/。 - Jason S

样本标准差，可以表示为 sqrt(S/(N-1))。 - Gregor y

同时，“S”现在也面临着“M”曾经遇到的增长问题。 - Gregor y

你可以去找高德纳问一下，也许他会给你寄一张2.56美元的支票。 - Jason S

Answer 2

这里有一种使用整数计算平均值的方法，避免了舍入误差和大中间值:

sum = 0
rest = 0
for num in numbers:
  sum += num / N
  rest += num % N
  sum += rest / N
  rest = rest % N

return sum, rest

- sepp2k

---

这基本上使用了多精度（双字）算术。我认为有一种方法可以优化，以减少除法（/或％）操作的数量，但我无法立即想起来。 - Jason S

通常的做法是在一个函数/单个操作中计算X/N和X%N。这是因为底层算法基本相同。 - MSalters

不，我更指的是： sum += (num + rest) / N; rest = (num + rest) % N;但这样做可能容易溢出。 - Jason S

C语言没有单独的函数或运算符来实现除法和取模，但是大多数优秀的编译器会在你同时使用两者时自动生成相应的指令或调用。 - Stephen Canon

只是一个小改进：由于“rest”永远不可能大于“2*(N-1)”，所以您可以将第二个除法/取模操作替换为一个简单的“if(rest > N)”分支（当然，这取决于平台分支或除法更快，还是编译器足够聪明地看到它自己）。实际上，为了避免溢出（如果“N>MAXINT/2”），每次“rest”变为正数时都可以从“rest”中减去“N”( 并在循环后决定是否要将“sum”四舍五入，具体取决于“rest”)。 - chtz

@Stephen： "C语言没有一个单一的函数或运算符来公开div/mod" - 除非我们计算div，它确实可以做到这一点 - 尽管不能用于浮点数。 - Eric

Answer 3

如果大整数是个问题...这样可以吗？

a/N + b/N+.... n/N

我理解你是在寻找其他方法或最优方法吗？

- Svetlozar Angelov

---

为什么？！如果a，b等是整数，这将给出一个错误的答案。如果它们是浮点数，我不确定，但我的直觉是，你会得到比只执行求和然后除以更多的舍入误差。在任何情况下，计算时间都会大大增加，而收益却值得怀疑。 - Jason S

Answer 4

如果数组是浮点数据，即使是“简单”算法也会受到舍入误差的影响。有趣的是，在这种情况下，将计算分块成sqrt（N）个长度为sqrt（N）的总和实际上可以减少平均情况下的误差（即使执行相同数量的浮点舍入）。

对于整数数据，请注意您不需要一般的“大整数”，如果您的数组中元素数量小于40亿（很可能），则只需要一个大于数组数据类型的32位整数类型。在这个稍微大一点的类型上执行加法基本上始终比在类型本身上执行除法或取模更快。例如，在大多数32位系统上，64位加法速度比32位除法/取模速度更快。随着类型变得更大，这种效应只会变得更加夸张。

- Stephen Canon

Answer 5

如果你使用 float，你可能会避免大整数：

def simple_mean(array[N]):
    sum = 0.0 # <---
    for i = 1 to N
       sum += array[i]
    return sum / N

- Nick Dandoulakis

Answer 6

Kahan算法（根据维基百科）具有比成对求和更好的运行时性能-O(n)-和O(1)的误差增长：

function KahanSum(input)
    var sum = 0.0
    var c = 0.0                  // A running compensation for lost low-order bits.
    for i = 1 to input.length do
        var y = input[i] - c     // So far, so good: c is zero.
        var t = sum + y          // Alas, sum is big, y small, so low-order digits of y are lost.
        c = (t - sum) - y // (t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
        sum = t           // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
        // Next time around, the lost low part will be added to y in a fresh attempt.
    return sum

其想法是单独从主求和中对浮点数的低位求和并进行校正。

- Gábor Bakos

Answer 7

在Jason S的解决方案的基础上构建，以找到加权平均值并控制S的增长。

使用之前给出的M查找算法以及加权平均值和总体标准差的聚合公式：

• Avg = Avg(W*X) / Avg(W)
• StDev = sqrt(Avg(W*X*X) / Avg(W) - Avg*Avg)

重写代码以找到三个运行平均值，然后在最后进行聚合计算。

function GetPopulationStats{
<#
.SYNOPSIS
Calculate the average, variance, and standard deviation of a weighted data set

.DESCRIPTION
Uses the Knuth method for finding means adapted by Jason S, to calculate the 
three averages required by the agregate statistical formulas

.LINK
https://dev59.com/mXM_5IYBdhLWcg3wgzdl#1346890
#>
   param(
      [decimal[]]$x #Data Points
     ,[decimal[]]$w #Weights
   )
   
   $N = $x.Length
   
   [decimal]$AvgW   = 0
   [decimal]$AvgWX  = 0
   [decimal]$AvgWXX = 0

   for($i=0; $i -lt $N; $i++){
      $AvgW   += ($w[$i]               - $AvgW)   / ($i+1)
      $AvgWX  += ($w[$i]*$x[$i]        - $AvgWX)  / ($i+1)
      $AvgWXX += ($w[$i]*$x[$i]*$x[$i] - $AvgWXX) / ($i+1)
   }

   [ordered]@{
      N     = $N
      Avg   = ($avg = $AvgWX / $AvgW)
      Var   = ($var = ($AvgWXX / $AvgW) - ($Avg * $Avg))
      StDev = SquareRoot $var
   }
}

如果你的编程语言类似于Windows PowerShell，那么你可能需要更高精度的[math]::sqrt()函数。

function SquareRoot([decimal]$a){
<#
.SYNOPSIS
Find the square-root of $a

.DESCRIPTION
Uses the McDougall-Wotherspoon variant of the Newton-Raphson method to 
find the positive zero of:
   f(x)  = (x * x) - a
   f'(x) = 2 * x

.NOTES
ToDo: using a fitted polynomial would likely run faster
#>
   $BiCycleX = $PrevX = 0; 
   $x  = $a/2 # guess
   $y  = ($x * $x) - $a
   $xx = $x
   $m  = $x + $xx

   $del = $x - $PrevX
   if($del -lt 0){ $del = -$del }
   
   $i = 0
   while($del -gt 0 -and $x -ne $BiCycleX){
      $BiCycleX = $PrevX;
      $PrevX    = $x;
      $x  = $x - ($y / $m)
      $y  = ($x * $x) - $a
      $xx = $x - ($y / $m)
      $m  = $x + $xx

      $del = $x - $PrevX
      if($del -lt 0){ $del = -$del }
      
      if(++$i -ge 50){
         throw ("invariant sanity fail on step {0}:`r`n   x_(n-1) = {1}`r`n   x_n     = {2}`r`n   delta   = {3:0.#e0}" -f $i,$PrevX,$x,$del)
      }
   }
   
   ($x + $PrevX) / 2
}

然而，如果你不需要加权方案，那么只需让i为所有w[i] = 1即可轻松解决。

最后，在代码上进行快速的合理性检查也是很有必要的。

describe 'tool spot-check' {
   context 'testing root calcs' {
      $TestCases = @(
         @{Value = 0; Expected = 0}
         @{Value = 1; Expected = 1}
         @{Value = 4; Expected = 2}
         @{Value = 9; Expected = 3}
         @{Value = 2; Expected = [decimal]'1.4142135623730950488016887242'}
         @{Value = (1e14-1); Expected = [decimal]'9999999.99999995'}
      )
      It 'finds the square root of: <Value>' -TestCases $TestCases {
         param($Value,$Expected)
         SquareRoot $Value | should be $Expected
      }
   }
   
   context 'testing stat calcs' {
      It 'calculates the values for 1 to 1000' {
         $x = 1..1000
         $w = @(1) * 1000
         
         $results = GetPopulationStats $x $w
         $results.N     | should be 1000
         $results.Avg   | should be 500.5
         $results.Var   | should be 83333.25
         $results.StDev | should be ([decimal]'288.67499025720950043826670416')
      }

      It 'calculates the values for a test data set' {
         $x = @(33,119,37,90,50,94,32,147,86,28,50,80,145,131,121,90,140,170,214,70,124)
         $w = @(207,139,25,144,72,162,93,91,109,151,125,87,49,99,210,105,99,169,50,59,22)
         
         $results = GetPopulationStats $x $w
         $results.N     | should be 21
         $results.Avg   | should be ([decimal]'94.54433171592412880458756066')
         $results.Var   | should be ([decimal]'2202.659150711314347179152603')
         $results.StDev | should be ([decimal]'46.93249567955356821948311637')
      }
   }
}

- Gregor y