サンプルサイズが大きいとき、単回帰分析における最小二乗推定量は近似的に正規分布に従う

要約

　サンプル・サイズが大きいとき、最小二乗法(OLS)によって計算される標本回帰係数は、近似的に正規分布に従います。つまり、単回帰分析における最小二乗推定量は、漸近正規性を持ちます。ただし、外生性といった条件を満たす必要があります。

　例えば、次図のようなモデルから生まれるサンプルを用いて、単回帰分析をして、標本回帰係数を求めたとしましょう。

　上のモデルに従うサンプル・サイズ1万の標本から得られる標本回帰係数のヒストグラムは、正規分布に近似します。

全体像

（１）問題の構造

　問いとして「サンプル・サイズnが大きいとき、最小二乗法で求めた標本回帰係数は近似的に正規分布に従うのか？」を設定します。この問いは

①サンプル・サイズnが大きいとは？

②最小二乗法で求めた標本回帰係数とは？

③近似的に正規分布に従うとは？

④母集団に課せられた条件とは？

⑤以上を踏まえて、サンプル・サイズnが大きいとき、最小二乗法で求めた標本回帰係数は近似的に正規分布に従うのか？

と分解できます。

（２）前提の選択

　前提として①②③④を選択します。

（３）論点の選択

　論点として⑤を選択します。

（４）付録一覧

冗長さを避けて、可読性を上げるために、以下の内容は付録に回します。

「Rコード」

前提

①サンプル・サイズnが大きい

　サンプル・サイズが大きいの定義を、100や1万と仮置きします。ここは適当です。

②最小二乗法で求めた標本回帰係数

　最小二乗法は、観測されたデータとモデルからの予測の差の二乗和を最小化することで、パラメータを推定します。統計学や機械学習などの分野で広く使用される方法です。R言語のlm関数にても採用されています。

③近似的に正規分布に従う

　「近似的に正規分布に従う」とは、得られたサンプルのヒストグラムがだいたい正規分布に従っていることとざっくり定義します。

④母集団に課せられた条件

　以下の4条件が正しいことを仮定します。

$$データ（X,Y）は独立で同一の分布に従い、単回帰モデルに従う$$

$$標本が無作為抽出$$

$$説明変数Xの分散が0でない$$

$$外生性：E(U|X)=0$$

　

　また、以下の2条件が正しいことを仮定しません。

$$均一分散：Var(U_i)=\sigma^2$$

$$誤差項の正規性：U_i　〜　N(0,\sigma^2)$$

　つまり、誤差項はXの条件つき期待値が0であれば、どんな分布に従っていても構いません。

方法

　まず、④を満たす母集団を一つ考えます。

　次に、「その母集団から標本抽出し、回帰分析して、標本回帰係数を求める」を5000回繰り返します。

　得られた5000個の標本回帰係数のヒストグラムと、正規分布を比較します。

　もし標本回帰係数のヒストグラムがだいたい正規分布であれば「(ある特定の場合ではあるが)サンプル・サイズnが大きいとき、最小二乗法で求めた標本回帰係数は近似的に正規分布に従う」と結論できます。

結果

（１）母集団を作る

単回帰モデル

　次の単回帰モデルを考えます。

$$Y=2+2X+U$$

$$Y：目的変数、X：説明変数、U：誤差項$$

三峰型の分布に従う誤差項

　下図のように明らかに正規分布に従っていない、誤差項U1を考えます。誤差項の母平均は0で、母分散は0.52です。黒は、正規分布N(0, 0.52)です。↓

不均一分散する誤差項

　さらに、不均一分散の要素を加えます。 U1は、説明変数X=1のときの誤差項とします。そして、説明変数X=xのときの誤差項Uxは↓

$$誤差項の期待値　E（U_x）=0$$

$$誤差項の分散　Var（U_x）=x \times Var（U_1）$$

標本抽出の例（n=500）

　すると、サンプル・サイズ500の場合、次のような標本が生成されます。

（２）回帰係数の分布

　標本生成と回帰分析を5000回繰り返して、回帰係数の推定値の分布を調べましょう。

n=1万

　サンプル・サイズn＝1万の場合を考えてみます。結果が下図です。正規分布に近似できると言えましょう。

n=100

　サンプル・サイズn＝100の場合を考えてみます。結果が下図です。正規分布に近似できていると言えましょう。

n=5

　逆に、サンプル・サイズが小さい場合（n＝５）も考えてみます。結果が下図です。正規分布に近似できているとは言えません。

考察

（1）結論

　サンプル・サイズnが大きいとき、最小二乗法で求めた標本回帰係数は近似的に正規分布に従います。

　ただし、以下の条件が必要です。

$$データ（X,Y）は独立で同一の分布に従い、単回帰モデルに従う$$

$$標本が無作為抽出$$

$$説明変数Xの分散が0でない$$

$$外生性：E(U|X)=0$$

（2）妥当性評価

前提評価

　母集団に課している条件、課していない条件が明確な点が、Goodです。

　サンプル・サイズが大きいの定義が適当な点が、Badです。

方法評価

　乱数発生によるモンテカルロシミュレーションを行い、視覚的に結果を表現したのは、Goodです。

　母集団一つの場合についてしか検討してないので一般性がない点が、Badです。

結論評価

　最小二乗推定量が漸近正規性を持つことはよく知られた性質であるので、Goodです。

　単回帰分析のみ検討しており、重回帰分析については検討できていないのが、Badです。

（3）意義

　推定量が正規分布に従うのならば、t検定を用いることができます。

　その点で誤差項の正規性を仮定しなくても、仮説検定できる可能性を切り開いた「最小二乗推定量の漸近正規性」は非常にパワフルな性質になります。

付録：Rコード

　シミュレーションは、R言語を用いて行いました。R言語については「しまうまのRでデータ分析入門」をご覧ください。

　美しい図示のため、パッケージとしてggplot2を用いました。

#パッケージの呼び出し。未インストールならインストール！
library(ggplot2) 

　誤差項の母集団の生成と図示↓

u_size <- 10000　
u0 <- rep(NA, u_size)
u1 <- rep(NA, u_size)
　
#三峰型の誤差項の「もと」を生成
for(i in 1:u_size){
  u0[i] <- ifelse(i < 3000,rnorm(1,-10,1) , 
                 ifelse(i < 6000, rnorm(1,-2,1),rnorm(1,7,2)))
}
　
#u1を作成。誤差項の「もと」の平均を0に調整＝外生性を仮定
for(i in 1:u_size){
  u1[i] <- 0.1*(u0[i] - mean(u0))　
}
　

#ヒストグラムの図示
data_u <- data.frame(u1)
　
ggplot(data_u, aes( x = u1)) + 
  geom_histogram( bins=50,aes( y = ..density.. ),colour="white", fill = "pink2")+
  stat_function(fun = dnorm, args=list(mean=mean(data_u$u1), sd=sqrt(var(data_u$u1))))+
  xlab("誤差項U") + 
  ylab ("確率分布") +
  ggtitle ("誤差項Uの母集団　　黒は、ピンクと同じ母平均、母分散の正規分布") +
  theme_grey(base_family = "HiraKakuPro-W3")

　標本生成と図示↓

sample_size0 <- 500 #サンプル・サイズを500
x<- rnorm(sample_size0, 15,3)
u_x <- rep(NA, sample_size0)
　
#標本生成
for(i in 1:sample_size0){ 
  xx <- x[i]
  uu <- sample(u1, 1, replace = TRUE)
  u_x[i] <- uu*xx
}
y <- 2 + 2*x+ u_x
　
#散布図作成
data0 <- data.frame(x, y)　
ggplot(data0, aes(x = x,y=y)) +　
  geom_point(aes() ) +
  scale_x_continuous(limits = c(0, 30),breaks = seq(0,30,5))+
  scale_y_continuous(limits = c(0, 80),breaks = seq(0,80,10))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)+
  xlab("説明変数X") + 
  ylab ("目的変数Y") +
  ggtitle ("サンプル・サイズ500の標本") +
  theme_grey(base_family = "HiraKakuPro-W3")

　シミュレーションで分布を図示↓

sample_size <- 100　#サンプル・サイズ
　
trial_number <- 5000　#試行回数
beta_1 <- rep(NA, trial_number)　
x <- rnorm(n=sample_size, mean =15, sd=3)
u_x <- rep(NA, sample_size) 
　
#5000回標本生成と回帰分析を繰り返す
for(j in 1:trial_number){ 
  for(i in 1:sample_size){ 
    xx <- x[i]
    uu <- sample(u1, 1, replace = TRUE)
    u_x[i] <- uu*xx
  }
  y <- 2 + 2*x+ u_x
  answer <- lm(y~x)
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
}
data00 <- data.frame(beta_1)　
　
#推定値の分布
ggplot(data00, aes( x = beta_1)) + 
  geom_histogram( bins=50,aes( y = ..density.. ),colour="white", fill = "pink2")+
  stat_function(fun = dnorm, args=list(mean=mean(data00$beta_1), sd=sqrt(var(data00$beta_1))))+
  xlab("回帰係数の推定値") + 
  ylab ("確率密度") +
  ggtitle ("最小二乗推定量の分布（n=100）　　黒は正規分布") +
  theme_grey(base_family = "HiraKakuPro-W3")