シミュレーションの方法

シミュレーションを行うために、データ生成を行う。私たちは真の関係を知っているので、その関係を利用する。

まず、標本サイズ \(n\) を決める。試しに、標本サイズを5にしてみよう。

n <- 5

次に、\(X\) の値を決める。とりあえず、\([-5, 5]\) からの一様分布から\(x\) ランダムに作ってみよう。

x <- runif(n, min = -5, max = 5)

続いて、\(y\) を生成する。真の関係は、\(Y = 2 + 0.8X + U\) である。

beta0 <- 2
beta1 <- 0.8

まず、\(U\) を作る。 \(U \sim \mbox{N}(0, \sigma^2)\) なので、誤差項の分散 \(\sigma^2\)（あるいは標準偏差 \(\sigma\)）を決める必要がある。ここでは、\(\sigma = 2\) としてみよう。

sigma <- 2

この標準偏差を使って、\(U\) をランダムに生成する。

U <- rnorm(n, mean = 0, sd = sigma)

これで、\(Y\) が生成できる。

y <- beta0 + beta1 * x + U

\(X\) の観測値 \(x\) と\(Y\) の観測値 \(y\) が手に入ったので、回帰分析を実行してみよう。

fit <- lm(y ~ x)
display(fit)

## lm(formula = y ~ x)
##             coef.est coef.se
## (Intercept) 3.40     0.64   
## x           0.72     0.20   
## ---
## n = 5, k = 2
## residual sd = 1.23, R-Squared = 0.81

\(\beta_0\)と\(\beta_1\)の推定値をそれぞれ \(b_0\)、\(b_1\) とすると、\(b_0\)=3.4009985、\(b_1\)=0.7246988 である。推定はどれくらい正確だろうか？

複数回のシミュレーション

1回だけのシミュレーションでは、その結果が偶然得られたものなのか、必然的な結果なのかの区別がつかない。そこで、上のシミュレーションを繰り返す。複数のコマンドを何度も実行するのは面倒なので、シミュレーション用の関数を function() で定義してしまおう。

単回帰 (simple regresion) のシミュレーションを行うので、関数名をsimple_reg にする。

simple_reg <- function(n, beta0 = 0, beta1 = 1, sigma = 1) {
  ## 単回帰のシミュレーションを実行するための関数
  ## 引数：n = 標本サイズ
  ##       beta0 = 真のy切片（規定値は0）
  ##       beta1 = 真の傾き（規定値は1）
  ##       sigma = 誤差項の標準偏差（規定値は1）
  
  # x を一様分布 U(-5, 5) から作る
  x <- runif(n, min = -5, max = 5)
  
  # u を正規分布 N(0, sigma^2) から作る
  u <- rnorm(n, mean = 0, sd = sigma)
  
  # 真のモデルからyを作る
  y <- beta0 + beta1 * x + u
  
  # 回帰分析を実行する
  fit <- lm(y ~ x)
  
  # beta の推定値を関数の出力として返す
  return(coef(fit))
}

function() で定義した関数は、return() で指定された対象を返し（これを戻り値または返り値, [return value] と呼ぶ）、終了する。ここで定義した関数は、切片と傾きの推定値を返す。

試しに、この関数を使ってみよう。私たちが実行したいのは、\(n=5\)、\(\beta_0 = 2\)、\(\beta_1 = 0.8\)、\(\sigma = 2\) の場合なので、次のようにする。

simple_reg(n = 5, beta0 = 2, beta1 = 0.8, sigma = 2)

## (Intercept)           x 
##    1.290141    0.486504

もう1度やってみよう。

simple_reg(n = 5, beta0 = 2, beta1 = 0.8, sigma = 2)

## (Intercept)           x 
##   2.3624474   0.8674906

もう1度やってみよう。

simple_reg(n = 5, beta0 = 2, beta1 = 0.8, sigma = 2)

## (Intercept)           x 
##   3.3825499   0.4624411

このように、実行する度に異なる結果が得られる。

これを繰り返し実行すれば、最小二乗法がどれくらい正確に推定を行えるか理解することができるはずである。 しかし、得られた結果が偶然の結果ではないと信じるためには、繰り返し回数を多くする必要がある。 たとえば、1,000回のシミュレーションを行う場合、上のように毎回コマンドを実行するのは面倒である。

そこで、forループを使ってシミュレーションを自動化しよう。まず、シミュレーション回数 n_sims を決める。

n_sims <- 1000

次に、シミュレーション結果を保存するための容器を用意する。私たちのシミュレーションでは、シミュレーションの繰り返し回数が n_sims 回、推定する母数（パラメタ, parameters）の数が2つなので、n_sims行、2列の行列を用意しよう。行列は、matrix() で作れる。 matrix() では、行列の要素と、行数 (nrow)、列数 (ncol) を指定する。ここでは空の容器を作りたいので、要素をすべて NA（欠測値）にした行列を作る。

result <- matrix(NA, nrow = n_sims, ncol = 2)

行列の最初の5行を確認してみよう。

result[1:5, ]

##      [,1] [,2]
## [1,]   NA   NA
## [2,]   NA   NA
## [3,]   NA   NA
## [4,]   NA   NA
## [5,]   NA   NA

要素がすべて NA になっていることがわかる。

わかりやすいように、行列の列に名前をつけておこう。

colnames(result) <- c("b0", "b1")

準備ができたので、forループ でシミュレーションを実行する。上で作った行列 result の \(i\) 行目に、\(i\)番目のシミュレーションの結果を保存する。

for (i in 1:n_sims) {
  result[i, ] <- simple_reg(n = 5, beta0 = 2, beta1 = 0.8, sigma = 2)
}

結果の最初の5行を確認してみよう。

result[1:5, ]

##            b0        b1
## [1,] 1.789458 0.1950699
## [2,] 1.925794 0.8248761
## [3,] 1.409744 1.6661211
## [4,] 2.777641 0.6866994
## [5,] 2.904388 0.6487705

シミュレーションの実行結果が保存されていることがわかる。

シミュレーションの結果を確認してみよう。 私たちが知りたいのは、回帰分析で、\(\beta_0=2\)、\(\beta_1 = 0.8\) がどれだけ正確に推定できるかということである。

まず、\(\beta_0\) の推定値である、\(b_0\) (b0) をヒストグラムにしてみよう。

res_data <- as.data.frame(result)  # 行列をデータフレームに変換する
hist_b0 <- ggplot(data = res_data, aes(x = b0)) +
  geom_histogram(color = "black") +
  labs(x = expression(b[0]), y = "度数") +
  geom_vline(xintercept = 2, color = "red")  # beta0 の真の値を示す
print(hist_b0)

ヒストグラムに加えられた赤い直線が真の値を示している。データ生成と推定を繰り返すと、推定がうまくいくこともあれば、そうでないこともあるということがわかる。分布の形に注目すると、分布の中心は真の値付近にあり、平均すると推定がうまくいっているように見える。

実際、b0の平均値は、

mean(res_data$b0)

## [1] 2.000733

であり、真の値である2に近い。

同様に、\(\beta_1\) の推定値である、\(b_1\) (b1) をヒストグラムにしてみよう。

hist_b1 <- ggplot(data = res_data, aes(x = b1)) +
  geom_histogram(color = "black") +
  labs(x = expression(b[1]), y = "度数") +
  geom_vline(xintercept = 0.8, color = "red")  # beta1 の真の値を示す
print(hist_b1)

ヒストグラムに加えられた赤い直線が真の値を示している。やはり、データ生成と推定を繰り返すと、推定がうまくいくこともあれば、そうでないこともあるということがわかる。分布の形に注目すると、分布の中心は真の値付近にあり、平均すると推定がうまくいっているように見える。

実際、b1の平均値は、

mean(res_data$b1)

## [1] 0.798378

であり、真の値である0.8に近い。

最後に、得られた回帰直線を図示してみよう。

plt <- ggplot(NULL) +
  geom_abline(intercept = res_data$b0, slope = res_data$b1, color = "gray") +
  geom_abline(intercept = 2, slope = 0.8, color = "royalblue") +
  xlim(-5, 5) + ylim(-3, 7)
print(plt)

この図は、推定された回帰直線をグレーで、真の回帰直線を青で描いている。回帰直線1つをランダムに選ぶと、その線は必ずしも真の関係を正しく捉えていない。しかし、平均してみると、真の直線の周りに推定された線が集まっていることがわかる。つまり、平均的には、回帰分析はうまくいきそうである。