３つの標準正規乱数の発生方法

嫌いな女性のタイプは靴の踵を踏んじゃってる人.

あっ, どーも僕です.

3つ手法による標準正規分布の発生

　こないだの本に従って, 一様乱数が得られるという仮定のもと, 3つ方法で標準正規乱数を発生させてみました.
　めんどくさいので, 理論も数式も書きませんが, Rのコードとグラフだけ紹介します. グラフに書かれている青い線はRのdnormをcurve関数で描いたものです. meanとsdは発生させた乱数の最尤推定量です.

中心極限定理

r_st_norm <- function (N = 100)
{
  X <- matrix (runif (N * 12), nrow = 12, ncol = N)
  Z <- colSums(X) - 6
  return(Z)
}

Box-Muller

r_stBM_norm <- function (N = 100)
{
  U1 <- runif (N)
  U2 <- runif (N)
  l_U1 <- sqrt (-2 * log (U1))
  Z1 <- l_U1 * cos (2 * pi * U2)
  Z2 <- l_U1 * sin (2 * pi * U2)
  return (list (z1 = Z1, z2 = Z2))
}

極座標法

r_stPC_norm <- function (N = 100)
{
  v1 <- NULL
  v2 <- NULL
  w <- NULL
  n <- 0
  while (n < N) {
    tv1 <- c(2 * runif (N * 2) - 1)
    tv2 <- c(2 * runif (N * 2) - 1)
    V <- tv1 ** 2 + tv2 ** 2
    index <- which (V < 1)
    v1 <- c(v1, tv1[index])
    v2 <- c(v2, tv2[index])
    w <- c(w, V[index])
    n <- length (w)
  }
  v1 <- v1[seq_len(N)]
  v2 <- v2[seq_len(N)]
  w <- w[seq_len(N)]
  c <- sqrt (- 2 * log (w) / w)
  z1 <- c*v1
  z2 <- c*v2
  return ( list (z1 = z1, z2 = z2))
}