単回帰分析における最小二乗推定量は、漸近正規性をもつのか?

 サンプル・サイズが大きいとき、ゆるい条件で、最小二乗推定量が正規分布に従うことを説明します。

 シミュレーションによる図解もあります。

要約

(1)問い

 サンプル・サイズnが大きいとき、最小二乗法で求めた標本回帰係数は正規分布に従うのでしょうか?

(2)結論

 外生性、大標本(n→無限大)などを仮定すると、最小二乗法で求めた標本回帰係数は、正規分布に従います。

$$\widehat{\beta_1}  〜   正規分布N$$

(3)意義

 以上は、一定の仮定のもとで、単回帰分析にて仮説検定ができることを意味します。

前提

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

 これに際して、次の定義・仮定・事実を用います。

(1)定義

 このモデルで用いる単回帰モデルは

$$Y_i=\beta_0+\beta_1X_i+U_i$$

です。概要は「単回帰モデル」をご覧ください。Yは目的変数、Xは説明変数、Uは誤差項、iはデータ番号です。

 この記事では標本回帰係数を、最小二乗法で求めた最小二乗推定量と考えます。

$$\widehat{\beta_1}:標本回帰係数$$

$$\beta_1:母回帰係数$$

(2)仮定する条件

 次の3つを仮定し

  1. データ(X,Y)は独立で同一の分布に従い、単回帰モデルに従う
  2. 標本が無作為抽出
  3. 説明変数Xの分散が0でない

 さらに

$$4.外生性:E(U_i|X_i)=0$$

を仮定します。

(3)仮定しない条件

 この記事では、以下を「仮定しません」。

$$系列相関なし:誤差項同士の共分散=0$$

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

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

(4)事実

最小二乗推定量

 仮定1〜3が満たされ、説明変数Xが確定しているとき、単回帰分析の最小二乗推定量は

$$\widehat{\beta_1}=\beta_1+ \frac{\sum\limits_{i=1}^n(x_i-\overline{x})U_i } {\sum\limits_{i=1}^n(x_i-\overline{x})^2 }$$

となることが知られています。詳しくは単回帰分析/ 最小二乗推定量とパラメーターの関係をご覧ください。

中心極限定理

 n個の確率変数Zが独立同一分布に従うとき、nが十分に大きいならば、確率変数Zの総和は特定の正規分布に従うことが知られています。詳しくは、中心極限定理をご覧ください。

方法

  • 計算する方法
  • シミュレーションする方法

の両方を行います。それぞれ「結果 <計算>」「結果 <シミュレーション>」をご覧ください。

結果 <計算>

 さて、

$$Z_i=\frac{(x_i-\overline{x})U_i } {\sum\limits_{j=1}^n(x_j-\overline{x})^2 }$$

とおきます。すると

$$標本回帰係数\widehat{\beta_1}=\beta_1+\sum\limits_{i=1}^n Z_i$$

です。

 さて、外生性の仮定より、UとXは独立なので、Zは独立同一分布に従います。nは十分に大きいので、中心極限定理により

$$\sum\limits_{i=1}^n Z_i 〜 正規分布N $$

です。したがって

$$\widehat{\beta_1}=\beta_1+\sum\limits_{i=1}^n Z_i 〜 正規分布N$$

です。つまり、最小二乗法で求めた標本回帰係数は、正規分布に従います。

結果 <シミュレーション>

(1)入力

 下図のように明らかに正規分布に従っていない、誤差項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)$$

 さて、この誤差項Uxを用いて標本(X, Y)を生成します。

$$Y_i = 2+ 2X_i + U_x$$

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

(2)出力

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

n=1万

 この記事で考えたのは、nが無限大の「漸近理論」「大標本理論」です。

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

n=100

 無限大と比べて、小さすぎる1万を用いたのに、最小二乗推定量は正規分布に近似できました。つまり、「大標本って具体的にn=いくら?」問題が発生しました。

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

n=10

 無限大と比べて、小さすぎる100を用いたのに、最小二乗推定量は正規分布に近似できました。つまり、「大標本って意外と小さい?」問題が発生しました。

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

n=3

 あまりによくできすぎていて、このシミュレーション手法自体の信頼性が疑問です。

 そこで、サンプル・サイズn=3の場合を考えてみます。結果が下図です。正規分布に近似できているとは言えません。しっかり「近似できない場合」のシミュレーションもできました。

考察

(1)結論

  • データ(X,Y)は独立で同一の分布に従い、単回帰モデルに従う
  • 標本が無作為抽出
  • 説明変数Xの分散が0でない
  • 外生性
  • サンプル・サイズが無限大

を仮定すると、最小二乗法で求めた標本回帰係数推は、漸近的に正規分布と従います。

 これは中心極限定理に帰着させた計算やシミュレーションによってわかりました。

(2)前提評価

  • 誤差項が系列相関していない
  • 誤差項が均一分散している
  • 誤差項が正規分布に従う

を仮定しなくて良いので、今回の仮定はゆるいと言えるでしょう。たしかに

  • サンプル・サイズが無限大

という大きな仮定を課していますが、シミュレーションではn=10程度でも正規分布にかなり近似できました。したがって、実際のところ、あまり気にしなくていい仮定と言えそうです。

(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)
  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)
    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")

コメント欄 お気軽にコメントをお寄せください!

タイトルとURLをコピーしました