どんな分布（注）からでもサンプリングされた標本の平均をさらに繰り返して標本の分布を考えたとき、その分布は正規分布になる。
というのでいろいろ分布を変えて試してみる。
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])) }
![小児科診療 2020年 07 月号 [雑誌]](https://m.media-amazon.com/images/I/51Tr9qH58sL._SL500_.jpg "小児科診療 2020年 07 月号 [雑誌]")
![整形外科 2020年 05 月号 [雑誌]](https://m.media-amazon.com/images/I/41d+MR7W5HL._SL500_.jpg "整形外科 2020年 05 月号 [雑誌]")
![臨床放射線 2019年 11 月号 [雑誌]](https://m.media-amazon.com/images/I/41YUVbhoLBL._SL500_.jpg "臨床放射線 2019年 11 月号 [雑誌]")
