ゼミの統計勉強会の資料を作っている時にたまたま出くわした話です。 今まで違いが分かっていなかったことが分かった、 という意味で若干テンションが上がり気味です。

問題設定

まず、ここに1から10まで0.1刻みで増えていくatomic vector Xが あったとします。

# X_1=1, X_2=1.1, X_3=1.2, ...
X <- seq(1, 10, 0.1)

この atomic vector に $y_i = \beta_0 + \beta_1x_i + \epsilon_i $ のような関数を当ててみます。 それぞれの index は個体を示していて、$x_1$ は 1.0 で$x_2$ は1.1、 と0.1ずつ増えていきます。そういう $x_i$ と $y_i$ の関係を示す式です。 なお、$ \epsilon_i $ は個体ごとの誤差で平均は0、標準偏差は1の標準正規分布を仮定しています。 また、切片 $\beta_0$ は0, 傾き $\beta_1$ は1です。

# function y
f.y <- function(x){
    b_0 <- 0
    b_1 <- 1
    e <- rnorm(1, mean=0, sd=1)
    y <- b_0 + b_1 * x + e
    return(y)
}

この時、Xに対して関数f.yを当てる方法は2種類あります。

Y <- sapply(X, f.y)
Y_b <- f.y(X)

これの挙動がどう変わるのか、というのが今回のミソです。

実際の挙動とグラフ

まず sapply の挙動ですが、 これはplotで確認できます。

plot(X, Y)

同様に Y_b <- f.y(X) の挙動も確認してみます。

plot(X, Y_b)

全然違いますね。何ででしょう。

理由

f(V) も sapply(V, f) も V の全ての要素を f() の引数に与えることに差はありません。 ただ、f(V) では e <- rnorm(1, mean=0, sd=1) が一度しか評価されないのに対し、 sapply(V,f) では e <- rnorm(1, mean=0, sd=1) が f() を V に適用するごとに評価されています。 その結果、誤差が1つかそれぞれにあるかが異なったということですね。

予想外の挙動なのでもっと仕様を知らないとならないなぁと思いつつ、 こういうのには今後も気をつけて行きたいです。