axjack's blog

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.

カイ二乗分布をRで図示する

ソース

# df個のN(0,1)な確率変数の二乗和 
f <- function(df) sum( rnorm(n=df,m=0,sd =1)^2 )

# 繰り返し回数
N <- 10000;

# カイ二乗分布をN回計算してヒストグラムを書く
showHistOfChiSqrd <- function (df){
cs <- sapply(rep(df,N),f)
hist(cs,
     freq = FALSE
     ,main=sprintf("Histogram of chi-squared distribution;\n degree of freedom = %d",df) 
     ,xlab = ""
     )
}

実行結果

df = 1 の場合

showHistOfChiSqrd(1) f:id:axjack:20181223194855p:plain

df = 3 の場合

showHistOfChiSqrd(3) f:id:axjack:20181223195048p:plain

df = 5 の場合

showHistOfChiSqrd(5) f:id:axjack:20181223195102p:plain

df = 10 の場合

showHistOfChiSqrd(10) f:id:axjack:20181223195115p:plain

実行結果の感想

  • *1教科書の図とほぼ同じとなった
  • 自由度が異なると分布の形が大きく異なることがわかった
  • Rのfunctionの書き方がわかった
  • sapplyは便利だ
    • for文使わずに済む
    • for文使わない方がめちゃくちゃ早い

ソースコード補足

sapply関数

cs <- sapply(rep(df,N),f)

 cs <- c( f(df), f(df), ... ,f(df) )

と同じ。

例:

> f1 <- function(x) 3*x + 2 ;
> y <- sapply(1:10, f1)
> y
 [1]  5  8 11 14 17 20 23 26 29 32

*1:統計検定2級対応 統計学基礎