標本の標準偏差の計算式

合計と二乗合計から標準偏差を計算したいことが時々あります。
頻繁には使わないのですが、なぜかぐぐってもこの式を見つけるのは苦労します。
自分でも計算できますが結構面倒です。
ということで、備忘録です。

定義

ここで、標本の標準偏差 u といっているのは次のことです。
（似たようなもので、n-1 ではなく n で割るものもあります。）

u = \sqrt{ \frac{1}{n - 1} \sum_{i=1}^{n}(x_i - \bar{x})^2 }

ここで、 $\bar{x}$ は、$x$ の平均。$\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i$

$u$ は、式を変形すると、次の様になります。

\begin{eqnarray}
u & = & \sqrt{ \frac{1}{n - 1} \sum_{i=1}^{n}(x_i - \bar{x} )^2 }
  \\
  & = & \sqrt{ \frac{1}{n - 1} \sum_{i=1}^{n}(x_i^2 - 2x_i\bar{x} + \bar{x}^2 ) }
  \\
  & = & \sqrt{ \frac{1}{n - 1} \sum_{i=1}^{n}
    \left\{
    x_i^2
    - 2x_i\frac{1}{n}\sum_{i=1}^{n}x_i
    + \left( \frac{1}{n}\sum_{i=1}^{n}x_i \right)^2
    \right\}
  }
  \\
  & = & \sqrt{ \frac{1}{n - 1}
    \left\{
    \sum_{i=1}^{n}x_i^2
    - 2\frac{1}{n}\sum_{i=1}^{n}x_i\sum_{i=1}^{n}x_i
    + n \left( \frac{1}{n}\sum_{i=1}^{n}x_i \right)^2
    \right\}
  }
  \\
  & = & \sqrt{
  \frac{n}{n - 1}
  \left\{
    \frac{1}{n} \sum_{i=1}^{n}x_i^2 
    -
    \left( \frac{1}{n} \sum_{i=1}^{n}x_i \right)^2
  \right\}
  }
  \\
  & = & \sqrt{
  \frac{n}{n - 1}
  \left\{
    \bar{ (x^2) } - (\bar{x})^2
  \right\}
  }
\end{eqnarray}

ここで、$\bar{x^2}$ は、$x^2$ の平均。$\bar{x^2} = \frac{1}{n}\sum_{i=1}^{n}x_i^2$

計算方法

サンプル数 n のデータ列で、合計を sum 二乗合計を sqsum とすると 平均と標本の標準偏差は、python では、次のように求めることができます。

# n : サンプル数 (n >= 2)
# sum : 合計
# sqsum : 二乗合計
# ave : 平均
# std : 標本の標準偏差
import math

ave = sum / n
std = math.sqrt((sqsum / n - ave * ave) * n / (n - 1))

注意点として、(sqsum / n - ave * ave) の引き算で有効桁数が失われる様な場合には誤差が大きくなってしまいます。

簡単な検算

#!/usr/bin/env python
import sys
import math

data = []
count = 100
if len(sys.argv) > 1:
    count = int(sys.argv[1])

# データの準備 -- data = [1, 2, ..., n-1]
for i in range(1, count+1):
    data.append(float(i))

#
# 変形した式での計算
#
sum = 0.0
sqsum = 0.0
for d in data:
    d = float(d)
    sum += d
    sqsum += d * d

n = float(len(data))
ave = sum / n
std1 = math.sqrt((sqsum / n - ave * ave) * n / (n - 1))

print("ave=%.3f" % ave)
print("method 1: std1=%.3f" % std1)

#
# 元の標本の標準偏差の定義の式での計算
#
sum2 = 0.0
for d in data:
    d = float(d)
    sum2 += (d - ave) * (d - ave)

std2 = math.sqrt(sum2 / (n - 1))
print("method 2: std2=%.3f" % std2)

結果は、以下のとおりで、method 1 と method 2 の値が同じでOKの様です。

$ ./std-test.py
ave=50.500
method 1: std1=29.011
method 2: std2=29.011