中心極限定理を確かめる

どんな分布（注）からでもサンプリングされた標本の平均をさらに繰り返して標本の分布を考えたとき、その分布は正規分布になる。
というのでいろいろ分布を変えて試してみる。
Rでは接頭語rdpqに続いて分布名があるので、標準で実装されている関数をapropos関数で取得すると、19個あるようだった。

 [1] "beta"     "binom"    "cauchy"   "chisq"    "exp"      "f"        "gamma"    "geom"     "hyper"   
[10] "lnorm"    "logis"    "nbinom"   "norm"     "pois"     "signrank" "t"        "unif"     "weibull" 
[19] "wilcox"  

19個の関数で必要なパラメータを適当に決め、30個サンプリングして平均値を求め、これを3000回繰り返すとする。
関数からサンプリングした標本の標本平均と標本分散はデータから決まるし、関数とパラメータから理論的な平均と分散も計算できるのでこれも計算する。
乱数発生関数に引数をいれるときに無理やりテキスト形式で

rnorm(n, mean, sd)

を実現するために、テキスト処理をしたあとにparse関数とeval関数で無理やり評価して関数を実行している。

ヒストグラムの設定が適当なせいなのか理論的な分布との重ね合わせが合ってないが、だいたいが平均値の分布は正規分布に従っている。
ここで、コーシー分布がおかしいことになっているが、コーシー分布は平均値が定義できない、となっているので、おかしいのが正しい。

library(stringr)
rdpq <- c("r", "d", "p", "q")
fs <- mapply(function(z) str_remove(apropos(sprintf("^%s", z)), sprintf("^%s", z)), rdpq)

f <- fs[[1]]
for(i in 1:(length(fs)-1)){
  for(j in (i+1):length(fs)){
    f <- intersect(f, intersect(fs[[i]], fs[[j]]))
  }
}

n <- 30 # ある分布からのサンプリング数
p <- ev <- replicate(length(f), vector("list"))
names(p) <- names(ev) <- f

p[[f[1]]] <- c(10, 20)    # beta
p[[f[2]]] <- c(20, 0.2)   # binom
p[[f[3]]] <- ""           # cauchy
p[[f[4]]] <- c(5)         # chisq
p[[f[5]]] <- c(2)         # exp
p[[f[6]]] <- c(4, 6)      # f
p[[f[7]]] <- c(0.5, 2)    # gamma
p[[f[8]]] <- c(0.3)       # geom
p[[f[9]]] <- c(5, 10, 12) # hyper
p[[f[10]]] <- c(0.3, 1.5) # lnorm
p[[f[11]]] <- c(10, 3.3)  # logis
p[[f[12]]] <- c(10, 0.25) # nbinom
p[[f[13]]] <- c(6, 3.5)   # norm
p[[f[14]]] <- c(5)        # pois
p[[f[15]]] <- c(10)       # signrank
p[[f[16]]] <- c(4)        # t
p[[f[17]]] <- c(2.8, 9.4) # unif
p[[f[18]]] <- c(4, 3)     # weibull
p[[f[19]]] <- c(4, 7)     # wilcox

niter <- 3000 # n回サンプリングするシミュレーションの試行回数
res <- replicate(niter,
colMeans(
mapply(function(z0, z1){
         eval(parse(text=sprintf("r%s(%d%s)", z0, n, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=",")))))
        }, f, p)
)
)

xq <- seq(0, 1, length=10000)
dens <- mapply(function(z0, z1){
         txt <- sprintf("q%s(%f%s)", z0, xq, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=",")))
         y <- mapply(function(w) eval(parse(text=w)), txt)
         d <- mapply(function(w) eval(parse(text=sprintf("d%s(%f%s)", z0, w, paste0(ifelse(all(z1==""), "", ","), paste0(z1, collapse=","))))), y)
         return(list(q=y, d=d))
        }, f, p, SIMPLIFY=FALSE)

ev[["beta"]] <- c(
  E=p[["beta"]][1]/(p[["beta"]][1]+p[["beta"]][2]),
  V=p[["beta"]][1]*p[["beta"]][2]/((p[["beta"]][1]+p[["beta"]][2])^2*(p[["beta"]][1]+p[["beta"]][2]+1))
  )

ev[["binom"]] <- c(
  E=prod(p[["binom"]]),
  V=p[["binom"]][1]*p[["binom"]][2]*(1-p[["binom"]][2])
  )

ev[["cauchy"]] <- c(
  E=NA,
  V=NA
  )

ev[["chisq"]] <- c(
  E=p[["chisq"]][1],
  V=2*p[["chisq"]][1]
  )

ev[["exp"]] <- c(
  E=1/p[["exp"]][1],
  V=(1/p[["exp"]][1])^2
  )

ev[["f"]] <- c(
  E=ifelse(p[["f"]][2] > 2, p[["f"]][2]/(p[["f"]][2]-2), NA),
  V=ifelse(p[["f"]][2] > 4, 2*p[["f"]][2]^2*(sum(p[["f"]])-2)/p[["f"]][1]/(p[["f"]][2]-2)^2/(p[["f"]][2]-4), NA)
  )

ev[["gamma"]] <- c(
  E=prod(p[["gamma"]]),
  V=p[["gamma"]][1]*p[["gamma"]][2]^2
  )

ev[["geom"]] <- c(
  E=1/p[["gamma"]][1],
  V=(1-p[["gamma"]][1])/(p[["gamma"]][1]^2)
  )

ev[["hyper"]] <- c(
  E=p[["hyper"]][3]*p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2]),
  V=p[["hyper"]][3]*p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2])*(1-p[["hyper"]][1]/(p[["hyper"]][1]+p[["hyper"]][2]))*(p[["hyper"]][1]+p[["hyper"]][2]-p[["hyper"]][3])/(p[["hyper"]][1]+p[["hyper"]][2]-1)
  )

ev[["lnorm"]] <- c(
  E=exp(p[["lnorm"]][1]+p[["lnorm"]][2]/2),
  V=exp(2*p[["lnorm"]][1]+p[["lnorm"]][2])*(exp(p[["lnorm"]][2])-1)
  )

ev[["logis"]] <- c(
  E=p[["logis"]][1],
  V=pi^2*p[["logis"]][2]^2/3
  )

ev[["nbinom"]] <- c(
  E=p[["nbinom"]][1]*(1-p[["nbinom"]][2])/p[["nbinom"]][2],
  V=p[["nbinom"]][1]*(1-p[["nbinom"]][2])/p[["nbinom"]][2]^2
  )

ev[["norm"]] <- c(
  E=p[["norm"]][1],
  V=p[["norm"]][2]
  )

ev[["pois"]] <- c(
  E=p[["pois"]][1],
  V=p[["pois"]][1]
  )

ev[["signrank"]] <- c(
  E=p[["signrank"]][1]*(p[["signrank"]][1]+1)/4,
  V=p[["signrank"]][1]*(p[["signrank"]][1]+1)*(2*p[["signrank"]][1]+1)/24
  )

ev[["t"]] <- c(
  E=ifelse(p[["t"]][1] > 1, 0, NA),
  V=ifelse(p[["t"]][1] > 2, 1/(1-2/p[["t"]][1]), NA)
  )

ev[["unif"]] <- c(
  E=mean(p[["unif"]]),
  V=diff(p[["unif"]])^2/12
  )

ev[["weibull"]] <- c(
  E=p[["weibull"]][2] * gamma(1 + 1/p[["weibull"]][1]),
  V=p[["weibull"]][2]^2 * (gamma(1 + 2/p[["weibull"]][1]) - (gamma(1 + 1/p[["weibull"]][1]))^2)
  )

ev[["wilcox"]] <- c(
  E=prod(p[["wilcox"]])/2,
  V=prod(p[["wilcox"]])*(sum(p[["wilcox"]])+1)/12
  )
ev <- lapply(ev, "*", c(1, 1/n))

# 分布からの推定と理論値の確認
m <- rowMeans(res)
v <- apply(res, 1, var)
a <- cbind(do.call(rbind, ev), m=m, v=v)

cols <- c("distribution"="orange", "theorem"="blue")
par(mfrow=c(length(f), 2), cex.main=2)
for(i in seq(f)){
  ml <- list(m[i], ev[[i]]["E"])
  vl <- list(v[i], ev[[i]]["V"])
  de <- mapply(function(z1, z2) mapply(qnorm, xq, z1, z2), ml, vl, SIMPLIFY=FALSE)
  y <- mapply(function(z0, z1, z2) mapply(dnorm, z0, z1, z2), de, ml, vl, SIMPLIFY=FALSE)
  yl <- c(0, max(unlist(y), na.rm=TRUE))
  hist(res[i,], probability=TRUE, xlab="mean", ylab="density", main=f[i], ylim=yl, cex.main=2)
  for(j in seq(y)){
    lines(de[[j]], y[[j]], col=cols[j], lwd=5)
  }
  legend("topright", legend=names(cols), col=cols, pch=15, cex=2)
  plot(dens[[i]]$q, dens[[i]]$d, type="l", lwd=5, xlab="", ylab="density", main=sprintf("%s distribution", f[i]))
}