正規分布に従う乱数発生の方法

ほぼ以下を元ネタにしているのでこちらの方が詳細に書いてあると思う。
60. 確率分布と乱数

3 変量正規分布に従う乱数を生成する関数 r3norm() を定義する

> r3norm <- function(mu, A, n) {
>    U <- svd(A)$u;  V <- svd(A)$v;  D <- diag(sqrt(svd(A)$d))
>    B <- U %*% D %*% t(V) # 行列 A の平方根
>    w <- c()
>    for (i in 1:n) w <- append(w, list(mu + B%*%cbind(rnorm(3))))
>    return(w)
> }
> mu3 <- cbind(c(1,1,1))   # 平均ベクトル(縦ベクトル)
> A3 <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
> w3 <- r3norm(mu3, A3, 2000)
> head(w3)
[[1]]
         [,1]
[1,] 3.760969
[2,] 1.506681
[3,] 2.271915
---- 以下省略

では、4変量正規分布に従う乱数を生成する関数 r4norm() を定義すると

> r4norm <- function(mu, A, n) {
>    U <- svd(A)$u;  V <- svd(A)$v;  D <- diag(sqrt(svd(A)$d))
>    B <- U %*% D %*% t(V) # 行列 A の平方根
>    w <- c()
>    for (i in 1:n) w <- append(w, list(mu + B%*%cbind(rnorm(4))))
>    return(w)
> }
> mu4 <- cbind(c(1,1,1,1)) # 平均ベクトル(縦ベクトル)
> A4 <- array(c(2,1,1,1,1,2,1,1,1,1,2,1,1,1,1,2), dim=c(4,4))
> w4 <- r4norm(mu4, A4, 2000)
> head(w4)
[[1]]
           [,1]
[1,]  0.9252319
[2,] -0.4436436
[3,] -0.3861779
[4,] -0.1728883
---- 省略

この４乱数が本当に正規分布を描いているか確認してみる。

# 一つ目のデータだけ取り出す
> y<-c()
> x<-c(1:2000)
> for (i in 1:2000) y<-append(y, w4[[i]][1])
> head(y)
[1]  0.9252319 -0.6451168  3.6471124 -1.1930477  0.4600912  1.0727863

# Summaryデータから平均を取得
> summary(y)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.48500  0.03779  0.99330  0.98500  1.96500  5.68700
> quantile(y) # クォンタイル点
         0%         25%         50%         75%        100% 
-3.48457453  0.03778559  0.99330414  1.96540774  5.68672752
# 標準偏差を取得
> variance <- function(x) var(x)*(length(x)-1)/length(x)
> sqrt(variance(y))
[1] 1.389023

> # グラフ描画
> plot(y,dnorm(y, mean=0.9850483, sd=1.38937), type="n")
> curve(dnorm(x, mean=0.9850483, sd=1.38937), type="l",add=T)

確かにちゃんと正規分布が描かれています。納得。