StatsBeginner: 初学者の統計学習ノート

文系出身で工学部(工学研究科)の教員をしてます。

不偏分散の分母がn-1であることの直観的な理解

 同じタイトルで数ヶ月前にエントリを書いていたのですが、よくよく考えると全然意味のある説明になってなかったので、まとめ直します。
 不偏分散が
 \displaystyle s^2 = \frac{1}{n - 1} \sum_{i=1}^{n} (x_i - \bar{x})^2
であることの、証明というか導出は、統計学の教科書を見れば書いてあるのですが、「なんで1引くの?」というのを直観的に理解したいと思ったことがある人は多いと思います。「データのばらつきは、独立したランダム(確率的)な変動が生じる情報空間の次元(広さの一種)を意味する自由度に応じて評価すべきもので、パラメータを1個(ここでは母平均)推定するとデータの変動に決定論的な仮定(拘束)が持ち込まれて自由度が1下がるから」というような説明で、直観的に腑に落ちるという人はいいのですが、そういう人はたぶん少数派でしょう。
 それで整理の仕方を考えていて、ひとまず以下のようにまとめましたが、これで腑に落ちる人もあまりいないかも知れません(笑)


 母平均と母分散を設定してやると、それにしたがって正規乱数を生成するような、乱数発生器を考えます。もちろん、観測値はすべてiid(独立同分布)とします。
 分散について考える時、ふつうはこういう数直線上の分布をイメージすると思います。これは母平均2、母分散1に設定しているのですが、標本平均は少しズレていて2より少し大きくなっていますね。



 ここで少しデータの見方を変えて、乱数発生器から生成されたサンプル(データセット)を、「サンプルサイズNを次元数とする空間」上にマッピングすることにします。イメージしづらいかも知れませんが、データポイント1つが1本の座標軸に対応しているような、N次元の空間を設けてやると、N個のデータからなる「データセット」そのものを、その空間上の1点として表現できるわけです。
 クラスの生徒が40人いるなら、生徒1人1人が座標軸となって、それが40本直交している、40次元の超空間を考えます。軸には0から100までの目盛りが売ってあり、たとえば算数のテストを実施したら、そのテストの「40人分の成績」を、この40次元空間上の1点としてプロットできます。国語の成績、社会の成績なども、同様にプロットできますし、算数のテストを何回も実施したなら、それも全部プロットできる。
 40次元の空間を図示することはできないので、3人の生徒しかいないという前提で、3次元の空間に50回のテストの結果をプロットするとこんな感じになります。



 この場合、観測値が「独立にばらつく方向」が、サンプルサイズ分(つまり3本)だけ存在すると言え、この捉え方の下では、「分散」とはすなわち1方向あたりの平均的なばらつき具合ということになります。乱数実験であれば母平均が分かるので、1つ1つの座標軸にその母平均の印を付けておくと、プロットされた1つのデータセットが「それぞれの軸に関して母平均からどれぐらいズレているか」を測ったものが母分散に一致し、それを軸の数だけ足し合わせると、総二乗誤差になります。
 そして、3個の観測値がすべて等しくなる点を結んだ線(2次元平面の場合でいうy=xの線)をマクロ経済学にならって45度線と呼ぶことにすれば(対角線でもいいが)、45度線からの距離が、そのサンプルの二乗誤差を表すことになりますね。



 ばらつきの総和としての総二乗誤差(偏差平方和)は、ばらつきの方向の数に比例します。つまり「分散×サンプルサイズ」が、総二乗誤差の期待値となります。しかしこれはあくまで、「乱数発生器に設定した母平均からの二乗誤差」を考える場合の話ですね。
 一方、サンプリングのたびに得られる「そのサンプルの標本平均からの二乗誤差」を考えるとすると、その期待値は上述の値よりも必ず減ります。なぜなら、じつは標本平均というのは理論上、「そこから測った二乗誤差を最小にする点」だからです。ちなみに中央値は、二乗しない「偏差」(の絶対値の和)を最小にする点です。この、「標本平均から測ると二乗誤差が最小になる」という話は、上述の特殊な空間で考えるよりも、元の数直線で考えたほうがイメージしやすいと思います。


 では、具体的にどれだけ「減る」のかというと、ちょうど「方向1本分」だけ減ります。したがって、総二乗誤差を方向の数で割って「1方向あたりの二乗誤差」(すなわち分散)に直すという処理をする際に、割る数をnではなくn-1にしておかないと、期待値が本来の「1方向あたりの二乗誤差」(分散)に一致しなくなるわけです。


 減る量が「ちょうど1方向分」だけであるのはなぜなのか。標本平均というのは、「そこから測った誤差がn-1個分だけ分かれば、残り1個のデータは自動的に決まってしまうような点」であり、これはすなわち、「標本平均から測ると、ばらつきの方向が1本減る」ということに等しい。このことを指して「残差の自由度が1下がっている」と言っていることになります。
 実際には、観測そのものが1個減るわけではないので、変動の方向も減ってるわけではないのですが、全体としてちょうど1本分減るように、すべての観測値が平均に引き寄せられます。……いや、じつはよく考えるとこの言い方にも語弊があって、標本平均というのはそこからデータが生成される「母」(母平均)ではなく、「本当の母がどこにいるかはわからないので、サンプルが得られるたびに仮にそこに母がいることにしよう」と決めている点に過ぎないので、正確には、平均からの偏差がトータルで方向1本分減るように、平均のほうが観測値に近づいていくというほうが、イメージとしては正しいです。


 「n-1で割る理由」を直観的に理解したいのなら、ここまでの説明で十分ですが、もう少し続けます。
 上述の話は、元のN次元空間のまま考えて、標本平均から測った場合の総二乗誤差の期待値が、乱数発生器に設定した母平均から測った総二乗誤差の期待値に比べて、方向1本分だけ減っているということです。一方、何回もサンプリングを繰り返してそのたびに残差(偏差)のベクトルを記録していくとすると、もちろん残差ベクトル自体はN次元のベクトルになるのですが、平均が特定の値になるように拘束を与えると(もちろん実際の観測ではそんな拘束はできませんが)、そのN次元空間の中に収まるn-1次元の超平面上に分布するようになります。
 つまり、基底を取り直せばn-1次元のデータとして再解釈できるということで、これが、「残差の自由度が1下がる」ということの直観的なイメージに対応します。
 下の図は、3次元空間上に分布するはずのデータに、平均が特定の値になるような拘束を与えた場合のグラフです。




 実際の計測では、サンプリングのたびに標本平均自体も変わるので、同じ平面上に分布するようにはならないのですが、いわば1回1回のサンプリングで得られる計測値のセットが、「3次元で動くものを2次元に押し込める」のと同じ程度に、自由を失うということですね。(ここ、意味がわかりにくいと思いますが。)


 「自由度」って、統計の勉強をするときに「なんかそんな概念あるらしい」程度で済ませている人も多いと思うのですが(うちの学生はほとんどそうです)、「自由度あたりのばらつきを評価する」という形で登場する、けっこう重要な概念だということを理解しておく必要があります。「独立に変動する情報の個数」(自由度)が情報のばらつきを捉える上で重要だということは、考えれば分かるはずですけど、考えなくても分析の実務はできてしまうので放置されがちな気がします。


 作図のためのRのコードは以下のとおりです。

set.seed(1)

n <- 100  # 点の数
x <- rnorm(n, mean=2, sd=1)

col_pt <- rgb(0.1, 0.4, 0.9, alpha=0.40)  # 青で半透明

plot(x, rep(0, n),
     pch=16,
     col=col_pt,
     cex=1.0,  # 点の大きさ
     ylim=c(-0.1, 0.1),
     yaxt="n", ylab="",
     xlab="Value",
     main="Samples on a Number Line")

abline(h=0, col="gray")

# 平均
abline(v=mean(x), col="red", lwd=2)
library(plotly)
set.seed(1)

T <- 50          # テスト回数(点の数)
mu <- 60          # 平均点
sigma <- 10       # 点数のばらつき(標準偏差)

# 3人分の得点(各行が1回のテスト、各列が生徒)
X <- matrix(rnorm(T * 3, mean = mu, sd = sigma), ncol = 3)
colnames(X) <- c("Student1", "Student2", "Student3")

# 0~100に丸める
X <- pmin(pmax(X, 0), 100)

plot_ly() |>
  add_markers(
    x = X[,1], y = X[,2], z = X[,3],
    marker = list(
      size = 2,
      color = "rgba(30, 30, 200, 0.35)"  # 青の透明
    )
  ) |>
  layout(
    title = "Each test = one point in 3D (3 students)",
    scene = list(
      xaxis = list(title = "Student 1 score", range = c(0, 100)),
      yaxis = list(title = "Student 2 score", range = c(0, 100)),
      zaxis = list(title = "Student 3 score", range = c(0, 100)),
      aspectmode = "cube"
    )
  )
library(plotly)
set.seed(1)

T <- 50
mu <- 60
sigma <- 10

# データ生成
X <- matrix(rnorm(T * 3, mean = mu, sd = sigma), ncol = 3)
colnames(X) <- c("Student1", "Student2", "Student3")
X <- pmin(pmax(X, 0), 100)

# 対角線用データ(0〜100)
diag_vals <- seq(0, 100, length.out = 100)

plot_ly() |>
  add_markers(
    x = X[,1], y = X[,2], z = X[,3],
    marker = list(
      size = 2,
      color = "rgba(30, 30, 200, 0.35)"
    )
  ) |>
  add_trace(
    x = diag_vals,
    y = diag_vals,
    z = diag_vals,
    type = "scatter3d",
    mode = "lines",
    line = list(
      color = "black",
      width = 5
    ),
    name = "45-degree line (x=y=z)"
  ) |>
  layout(
    title = "Each test = one point in 3D (3 students)",
    scene = list(
      xaxis = list(title = "Student 1 score", range = c(0, 100)),
      yaxis = list(title = "Student 2 score", range = c(0, 100)),
      zaxis = list(title = "Student 3 score", range = c(0, 100)),
      aspectmode = "cube"
    )
  )
library(plotly)
set.seed(0)

N <- 1500
mu <- 2
sigma <- 1

X <- matrix(rnorm(N*3, mean=mu, sd=sigma), ncol=3)
colnames(X) <- c("X1","X2","X3")

R <- X - rowMeans(X)

grid <- seq(min(R), max(R), length.out=25)
z_mat <- outer(grid, grid, function(x,y) -(x+y))

# ふつうにプロット
p1 <- plot_ly() |>
  add_markers(x=X[,1], y=X[,2], z=X[,3],
              marker=list(size=1.5, opacity=0.35)) |>
  layout(title="Raw samples (R^3)",
         scene=list(
           xaxis=list(title="X1"),
           yaxis=list(title="X2"),
           zaxis=list(title="X3")
         ))

# 中心化した場合
p2 <- plot_ly() |>
  add_markers(x=R[,1], y=R[,2], z=R[,3],
              marker=list(size=1.5, opacity=0.5)) |>
  add_surface(x=grid, y=grid, z=z_mat,
              opacity=0.15, showscale=FALSE) |>
  layout(
    title="Centered residuals (R1+R2+R3=0)",
    scene=list(
      xaxis=list(title="R1"),
      yaxis=list(title="R2"),
      zaxis=list(title="R3"),
      aspectmode="cube"
    )
  )

p1  # 表示する
p2  # 表示する


(参考:このブログのおすすめ記事一覧はコチラ