正規分布からカイ二乗分布導く
よく知られる確率分布は互いに関連しており、ある分布からある分布を導出することができる。
カイ二乗分布は標準正規分布からランダムに得られた数値の2乗値の合計の分布である。
カイ二乗分布は自由度によって異なるが、ここでの自由度kは、標準正規分布の数と一致する。
すなわち、3つの標準正規分布からそれぞれn個の数値をランダムに取り出し、その二乗値の合計値の分布は、自由度3のカイ二乗分布となる。
なお、ここで合計というのは、例えば、3つの標準正規分布、A、B、Cがあり、それぞれランダムに1回（n=1）ずつ取り出した時、標準正規分布Aから取り出した数値の二乗をX^2A、標準正規分布Bから取り出した数値の二乗をX^2B、標準正規分布Cから取り出した数値の二乗をX^2Cとした時の、これらの合計のことである（=X^2A+X^2B+X^2C）。
よって、「3つの標準正規分布からn個の数値をランダムに取り出し、その二乗値の合計値」というのはn個の値が手元に得られることを意味する。
よって、自由度kのカイ二乗分布に従うn個の変数は、以下のように表現できる。なお、ここで、X^2は標準正規分布kからランダムに得られる値を2乗したものである。
実際にRで実験してみた。
なお、n=100,000で、自由度は1~9のカイ二乗分布を図示することを目的とした。
その結果とコードは以下の通り。
n<-100000 #n個
z<-matrix(nrow=n,ncol=9) #収納スペース
for(i in 1:9){
x<-rnorm(n) #標準正規分布からn個の乱数発生
x<-x^2 #n個の乱数をそれぞれ2乗
z[,i]<-x #収納
}
library(ggplot2)
g<-list() #収納
for(i in 1:9){
data<-data.frame(z)
data<-data[,1:i] #1~i列の乱数の2乗値を抽出
if(i == 1){ #以下のapply関数では列が1の時エラーとなるのでi=1の時は別にした
x2<-data
g[[i]]<-ggplot()
g[[i]]<-g[[i]]+geom_histogram(mapping=aes_string(x=x2),binwidth=0.1)
g[[i]]<-g[[i]]+xlim(0,30)+ylim(0,4000)
g[[i]]<-g[[i]]+labs(x="chi.square")
}
else{
x2<-apply(data,1,sum) ##1~i列の乱数の2乗値を合計
g[[i]]<-ggplot()
g[[i]]<-g[[i]]+geom_histogram(mapping=aes_string(x=x2),binwidth=0.1)
g[[i]]<-g[[i]]+xlim(0,30)+ylim(0,4000)
g[[i]]<-g[[i]]+labs(x="chi.square")
}
}
library(gridExtra)
g<-grid.arrange(g[[1]],g[[2]],g[[3]],g[[4]],g[[5]],g[[6]],g[[7]],g[[8]],g[[9]])
いい感じに図示できた。
同じ条件で直接カイ二乗分布から乱数を発生させて、それをヒストグラムにした場合も出してみる。
その結果とコードは以下の通り。
n<-100000
chi.x<-matrix(nrow=n,ncol=9)
for(i in 1:9){
chi.x[,i]<-rchisq(n=n,df=i) #自由度iのカイ二乗分布からn個の乱数発生させる
}
library(ggplot2)
g<-list()
for(i in 1:9){
data<-data.frame(chi.x)
x2<-data[,i]
g[[i]]<-ggplot()
g[[i]]<-g[[i]]+geom_histogram(mapping=aes_string(x=x2),binwidth=0.1)
g[[i]]<-g[[i]]+xlim(0,30)+ylim(0,4000)
g[[i]]<-g[[i]]+labs(x="chi.square")
}
library(gridExtra)
g<-grid.arrange(g[[1]],g[[2]],g[[3]],g[[4]],g[[5]],g[[6]],g[[7]],g[[8]],g[[9]])
予想通りほとんど同じ結果になった。
