確率密度関数の実験(1)カイ２乗分布

/Users/kazuo/Dropbox/RStudio/社会統計学/確率密度-練習１カイ２乗.Rmd

library(tidyverse)
library(showtext)
showtext_auto(TRUE)

1 概要

• 本レポートは、dnorm, pnorm, qnorm などの確率分布に関係するfunctionを理解するために、グラフを描き、教科書の巻末についている確率分布表を作成してみる。参考になるのは、(奥村 2016)、(Bruce et al. 2018)、(津島, 山口, and 田邊 2014)、(中澤港 2003)などである。functionの一覧的解説は、(船尾, n.d.)にある。

2 カイ2乗分布に関する基本的な関数

• 以下にグラフを表示するが、dchisq、pchisqの横軸はカイ二乗値である。そして、pchisqは、dchisq の累積値となっており、積分値である。
• qnormは、pchisqとは逆関数の関係になっており、ある分位点（積分値）になるカイ二乗値がいくつかを出力する。

2.1 確率密度関数（dchisq）を使う

カイ２乗分布に関連するfunctionをみていく。dchisqは、引数にx=カイ二乗値、df=自由度をとる。

x = c(1:2500)/250
p <- ggplot()
p <- p + geom_line(aes(x = x,y=dchisq(x,1)))
p <- p + ylim(c(0,1)) + ylab("dchisq(x,df)") + xlab("x = χ2値")
p <- p + ggtitle("自由度1のχ2分布の確率密度")
p

## Warning: Removed 34 row(s) containing missing values (geom_path).

2.2 参考：base graphicでのグラフ

x = c(1:1500)/250
y = dchisq(x,1) 
plot(x,y,type='s',main="df=1 のカイ2乗分布の確率密度関数",ylim = c(0,3))
abline(h=0)

2.3 累積確率密度関数（pchisq）を使う

\(\chi ^2\)値が、3.8415 となる時のα値はいくつになるか。おなじみの「0.05」になる値を求めたいののだが、カイ二乗値=0からの累積が面積なので、pchisq(3.8415,1) の値は、0.95になるので、それを1から引くと「0.05」になるという関係である。

1- pchisq(3.8415,1) # χ2値= 3.8415 の時の斜線部分

## [1] 0.04999877

x = c(1:2500)/250
p <- ggplot()
p <- p + geom_line(aes(x = x,y=pchisq(x,1)))
p <- p + ylim(c(0,1)) + ylab("pchisq(x,df)") + xlab("x = χ2値")
p <- p + ggtitle("自由度1のχ2分布の累積確率密度")
p

2.4 分位点関数 qchisq を使う

• 分位点(累積確率密度（つまり面積）がある値になるχ2値)を取得する。

確率密度関数のグラフで「ある値」から右側の部分の面積が「0.05」になるχ２値は、dchisq()を0からそのχ２値までの範囲の面積（積分）を１から引いたものである。 分位点関数でそれが取得できる。

qchisq(1 - 0.05,1) # 知りたいのは（おなじみの）0.05になるχ2値 求めるのは、1-0.05になるχ2値。

## [1] 3.841459

2.5 これを使えば「表」をつくれる

χ2乗分布表にあるような α 値になるχ2値を求めるには、以下のようにすればよい。自由度を１〜10で計算している。

N = 10
data.frame(df=1:N,
           A.1=qchisq(1-0.1,1:N),
           A.05=qchisq(1-0.05,1:N),
           A.025=qchisq(1-0.025,1:N),
           A.01=qchisq(1-0.01,1:N),
           A.005=qchisq(1-0.005,1:N),
           A.001=qchisq(1-0.001,1:N)
          ) %>% knitr::kable(digits = 4)

df	A.1	A.05	A.025	A.01	A.005	A.001
1	2.7055	3.8415	5.0239	6.6349	7.8794	10.8276
2	4.6052	5.9915	7.3778	9.2103	10.5966	13.8155
3	6.2514	7.8147	9.3484	11.3449	12.8382	16.2662
4	7.7794	9.4877	11.1433	13.2767	14.8603	18.4668
5	9.2364	11.0705	12.8325	15.0863	16.7496	20.5150
6	10.6446	12.5916	14.4494	16.8119	18.5476	22.4577
7	12.0170	14.0671	16.0128	18.4753	20.2777	24.3219
8	13.3616	15.5073	17.5345	20.0902	21.9550	26.1245
9	14.6837	16.9190	19.0228	21.6660	23.5894	27.8772
10	15.9872	18.3070	20.4832	23.2093	25.1882	29.5883

3 dchisqとpchisqの関係を理解する：確率密度関数(dchisq(x,df))を0〜ある値まで積分する。

• Rで積分するのは、integrate()
• これが、累積確率密度関数（pchisq）の値と同じになることを確認する。

fchi3 <- function(x){dchisq(x,3)} # 自由度3で0〜7.813の範囲を積分する
integrate(fchi3,0,7.8147)

## 0.9499994 with absolute error < 8.3e-09

# 以下のように書くと一発
integrate(function(x){dchisq(x,3)},0,7.8147)

## 0.9499994 with absolute error < 8.3e-09

http://www3.u-toyama.ac.jp/kkarato/2019/statistics/table/chisq.pdf

3.1 累積確率密度関数 pchisqを使って自由度1で、α=0.1, 0.05, 0.025, 0.01 となるχ2値を求める

3.1.1 力尽くで数値計算版

• χ2値を、小刻み（1/N）で変化させて、ターゲットのα値になるまでwhileでループさせている

get_chi2 <- function(df=1,alpha =0.05,N=10000){
  chi2 <- 1/N
  p <- 1
  while (p > alpha) {
    p <- 1-pchisq(chi2,df)
    chi2 <- chi2 + 1/N
  }
return(chi2)
}

chi2_table <- NULL
for(i in 1:10){
  chi2_table <- bind_rows(chi2_table,
                          c(i,a0.1=get_chi2(df=i,0.1),
                              a0.05=get_chi2(df=i,0.05),
                              a0.025=get_chi2(df=i,0.025),
                              a0.01=get_chi2(df=i,0.01)
                            ))
}
chi2_table %>% rename(df=...1)

## # A tibble: 10 x 5
##       df  a0.1 a0.05 a0.025 a0.01
##    <dbl> <dbl> <dbl>  <dbl> <dbl>
##  1     1  2.71  3.84   5.02  6.63
##  2     2  4.61  5.99   7.38  9.21
##  3     3  6.25  7.81   9.35 11.3 
##  4     4  7.78  9.49  11.1  13.3 
##  5     5  9.24 11.1   12.8  15.1 
##  6     6 10.6  12.6   14.4  16.8 
##  7     7 12.0  14.1   16.0  18.5 
##  8     8 13.4  15.5   17.5  20.1 
##  9     9 14.7  16.9   19.0  21.7 
## 10    10 16.0  18.3   20.5  23.2

3.2 累積確率密度関数（積分）を使って力づくでやった計算は、分位点関数::qchisq()で簡単に実現できる！

あっけなかった….。

alpha_list <- 1 - c(0.1,0.05,0.025,0.01,0.005,0.001) 
chisq_tbl <- NULL
for(df in 1:30){
  chisq_tbl <- bind_cols(chisq_tbl,qchisq(alpha_list,df))
} 
t(chisq_tbl) -> .tbl
colnames(.tbl) <- c(0.1,0.05,0.025,0.01,0.005,0.001) 
.tbl %>% as_tibble %>% mutate(df=1:30) %>% select(df,1:6) %>% knitr::kable(digits = 4)

df	0.1	0.05	0.025	0.01	0.005	0.001
1	2.7055	3.8415	5.0239	6.6349	7.8794	10.8276
2	4.6052	5.9915	7.3778	9.2103	10.5966	13.8155
3	6.2514	7.8147	9.3484	11.3449	12.8382	16.2662
4	7.7794	9.4877	11.1433	13.2767	14.8603	18.4668
5	9.2364	11.0705	12.8325	15.0863	16.7496	20.5150
6	10.6446	12.5916	14.4494	16.8119	18.5476	22.4577
7	12.0170	14.0671	16.0128	18.4753	20.2777	24.3219
8	13.3616	15.5073	17.5345	20.0902	21.9550	26.1245
9	14.6837	16.9190	19.0228	21.6660	23.5894	27.8772
10	15.9872	18.3070	20.4832	23.2093	25.1882	29.5883
11	17.2750	19.6751	21.9200	24.7250	26.7568	31.2641
12	18.5493	21.0261	23.3367	26.2170	28.2995	32.9095
13	19.8119	22.3620	24.7356	27.6882	29.8195	34.5282
14	21.0641	23.6848	26.1189	29.1412	31.3193	36.1233
15	22.3071	24.9958	27.4884	30.5779	32.8013	37.6973
16	23.5418	26.2962	28.8454	31.9999	34.2672	39.2524
17	24.7690	27.5871	30.1910	33.4087	35.7185	40.7902
18	25.9894	28.8693	31.5264	34.8053	37.1565	42.3124
19	27.2036	30.1435	32.8523	36.1909	38.5823	43.8202
20	28.4120	31.4104	34.1696	37.5662	39.9968	45.3147
21	29.6151	32.6706	35.4789	38.9322	41.4011	46.7970
22	30.8133	33.9244	36.7807	40.2894	42.7957	48.2679
23	32.0069	35.1725	38.0756	41.6384	44.1813	49.7282
24	33.1962	36.4150	39.3641	42.9798	45.5585	51.1786
25	34.3816	37.6525	40.6465	44.3141	46.9279	52.6197
26	35.5632	38.8851	41.9232	45.6417	48.2899	54.0520
27	36.7412	40.1133	43.1945	46.9629	49.6449	55.4760
28	37.9159	41.3371	44.4608	48.2782	50.9934	56.8923
29	39.0875	42.5570	45.7223	49.5879	52.3356	58.3012
30	40.2560	43.7730	46.9792	50.8922	53.6720	59.7031

3.3 qchisq を図示してみる

• pchisqのx、yが入れ替わっていることがわかる。

x = c(0:250)/250
p <- ggplot()
p <- p + geom_line(aes(x = x,y=qchisq(x,1)))
p <- p + ylab("qchisq(x,df)") + xlab("x = 累積値") # ylim(c(0,1))
p <- p + ggtitle("自由度1のχ2分布の分位点")
p

4 まとめ

4.1 確率密度関数（dchisq）、累積確率密度関数（pchisq）、そして分位点関数（qchisq）の関係

ラップ関数を用いて網掛け部分を描画する方法は (Chang 2019)p308-310を参照

dchisq_limit <- function(x){　# ラップ関数を用いて網掛け部分を描画する方法は p308-310を参照
  y <- dchisq(x,df=1)
  y[x < qchisq(0.9,1)] <- NA
  return(y)
}

p <- ggplot(data.frame(x=c(0,6)),aes(x = x))
p <- p + stat_function(fun = dchisq_limit, geom ="area",fill= "blue",alpha=0.2) + 
         stat_function(fun=dchisq,args = list(df=1)) + ylim(c(0,1)) +
         stat_function(fun=pchisq,args = list(df=1),colour="blue") 
p <- p + annotate("text", x = 1.3, y = 0.4, parse=FALSE,size=4,label="確率密度関数\ndchisq(x,df=1)") + 
  annotate("text", x = 0.7, y = 0.125, parse=TRUE,size=4,label="ここの面積0.9") + 
  annotate("text", x = qchisq(0.9,df=1), y = 0.2, parse=FALSE,size=4,label="分位点関数\nqchisq(0.9,df=1)\n=2.705543") +
  annotate("text", x = qchisq(0.9,df=1)+0.05, y = 0.1, parse=FALSE,size=4,label="↓") + 
  annotate("text", x = 4, y = 0.1, parse=FALSE,size=4,label="網掛けの面積は、1−0.9") +
  annotate("text", x = 3, y = 0.8, parse=FALSE,size=4,label="累積確率密度関数\npchisq(x, df=1)\n(確率密度関数の積分)") #+
#  annotate("text", x = 4.5, y = 0.75, parse=TRUE,size=4,label="integral(dchisq(x,df=1)*dx,0,chi)")
p <- p + ggtitle("カイ２乗分布（自由度１）で見る確率密度関数・累積確率密度関数・分位点関数")
p 

分位点関数（qchisq）と確率密度関数（dchisq）の関係は以下の通り。\[\int^{\chi^{2}}_0 dchisq(x,df=1)dx\]

4.2 自由度１〜4のカイ2乗分布の確率密度 ggplot2 で描画してみる

x = c(1:2500)/250
p <- ggplot() + geom_line(aes(x = x,y=dchisq(x,1)))
p <- p + geom_line(aes(x=x,y=dchisq(x,2)),colour="red")
p <- p + geom_line(aes(x=x,y=dchisq(x,3)),colour="blue")
p <- p + geom_line(aes(x=x,y=dchisq(x,4)),colour="green")

p <- p + ylim(c(0,1)) + ylab("dchisq(x,df)")
p <- p + ggtitle("自由度1~4のカイ二乗分布の確率密度")
p

## Warning: Removed 34 row(s) containing missing values (geom_path).

参考文献

Bruce, Peter C, 1953-, Andrew Bruce 1958-, 黒川利明(1948-), and 大橋真也. 2018. データサイエンスのための統計学入門: 予測、分類、統計モデリング、統計的機械学習とRプログラミング.

Chang, Winston. 2019. Rグラフィックスクックブック 第2版 ―ggplot2によるグラフ作成のレシピ集. Translated by 石井弓美子, 河内崇, and 瀬戸山雅人. 第2版 ed. 東京; 東京: オライリージャパン.

中澤港. 2003. Rによる統計解析の基礎. 東京: ピアソン・エデュケーション.

奥村晴彦. 2016. Rで楽しむ統計.

津島昌寛, 山口洋, and 田邊浩. 2014. 数学嫌いのための社会統計学.

船尾. n.d. “R-Tips.” http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html.