思いついたこと実験の記録
また気になる病にはまった。綺麗に抜けた感はない。。。
ABテストとピアソンのカイ二乗検定について知りたいだけの人は真ん中あたりで引き返してください
話す事
・ピアソンのカイ二乗値が期待値を分母に使っているのか？
・ピアソンのカイ二乗値がカイ二条分布に従うのが納得いかなかったのでシミュレーション
・ついでに検定ためす
・二項分布を変形することで、正規分布、さらに標準正規分布に変形できるか試す
・ピアソンのカイ二乗値が標準正規分布の二乗の合計になれば納得がいく(過去記事でやったので)
二項分布は正規分布に近似できる
・あたりかはずれか
・表か裏か
・感染するかしないか
など、二択の結果がある確率(割合)で起こるような試行「ベルヌーイ試行」
ベルヌーイ試行を複数回繰り返すことで考えられる分布「二項分布」
nが試行回数、pが表の出るパラメータ(確率)とする
\begin{align}
{}_n C_{x} p^x (1-p)^{n-x} \tag{1}
\end{align}
二項分布は試行回数が大きくなることで、真のパラメータと試行回数から計算できる値を中心とするような正規分布に近似できる
(スターリングの近似を使って、中心極限定理の特別な場合として二項分布は正規分布の関数に近似する-ド・モアブル-ラプラスの定理)
ここで思った
正規分布に近似できて、平均と分散が求まるならば標準化して「標準正規分布」にも変換できるし
標準正規分布にできれば二乗して合計することで「カイ二条分布」に従うのではないか？
というのも、ピアソンの定理を知らない状態でクロス集計表に対するピアソンの独立性検定を読んだ時に、
「ピアソンのカイ二乗値」は「カイ二乗値」とは別物である事を知り、
なぜカイ二乗検定でなく、ピアソンの統計量を使うのだろう？
と思ったので、確かめてみる
まず、「ピアソンのカイ二乗値」が自由度(α-1)×(β-1)のカイ二乗分布に従うことを確かめてみる
(しばらくコードを書かないと書き方も忘れる)
サイトAとBがあり、AとBの広告のクリック率の違いを確認するような場合を例にして確かめてみる。
こんなクロス集計表があったとする。
シミュレーションデータなので、クリック率は0.8でAもBも変化無いものとする。
n=1000
click = rbinom(2,n,0.8)
not_click = 1000-click
df = data.frame(click=click,not_click=not_click)
rownames(df) <- c("A","B")
df
click not_click
A 804 196
B 798 202
とりあえず独立性の検定のために「ピアソンのカイ二乗値」計算をする
rowSumsやcolSumsで表側や列合計を作る
以下が二つのサイトに違いが無いと考えたときの潜在的なクリック率(閲覧者数に対するクリック数の期待値 = 期待度数)になる
expected_frequency = colSums(df) * ( rowSums(df) / sum(rowSums(df)) )
nが十分に大きい二項分布では、正規分布に近似できる。
B(n,p) = N(np,np(1-p))
ここでピアソンさんは
「観測データからではバラツキってわからないよね。観測している範囲の値を使って近似できる分布を見つけられないかな？」
ということで「ピアソンのカイ二乗値」が生まれた。
\chi^2(df=(\alpha-1)*(\beta-1)) = \frac{(観測値 - 期待度数)^2}{期待度数} \\
= \frac{(X - np)^2}{np} \\
でも正規分布に近似できるといって、分散も示されているんだから
このまま分散に置き換えてやったら標準正規分布が作れて、カイ二条分布に近似できるんじゃない？
と思ったのが本記事の始まり。
\frac{(X - np)}{\sqrt{np(1-p)}} = N(0,1) ??? \\
z ～ N(0,1) \\
z_{1}^2 + z_{2}^2 + ・・・ +z_{k}^2 = \chi^2(df=k)
期待度数を観測値から引き算したら、期待度数から実際の観測値がどれだけ外れているかわかるし、
col_1 = (df[1] - expected_frequency[1])^2 / expected_frequency[1]
col_2 = (df[2] - expected_frequency[2])^2 / expected_frequency[2]
Pearson_chi <- sum(col_1[1,1] + col_1[2,1] + col_2[1,1] + col_2[2,1] )
Pearson_chi
0.1129242
ピアソンさんの話だとこれが自由度(行数-1)*(列数-1)の自由度のカイ二乗分布に従うってことらしいけど・・・
sum(col_1[1,1] + col_1[2,1] + col_2[1,1] + col_2[2,1] )
2行2列なので自由度は1
自由度1のカイ二乗分布では、カイ二乗値（検定統計量）が3.841となる時が有意水準5%の値(p-値0.95)である
(1%水準では「6.63」)
qchisq(0.95,df=1)#3.841
pchisq(3.841,1)#p値
curve(dchisq(x, 1),xlim=c(0,4),ylim=c(0,3))
abline(v=qchisq(0.95, (2-1)*(2-1)),col="red")
abline(v=Pearson_chi)
text(2, 2, labels=paste0("95% = ",round(qchisq(0.95, (2-1)*(2-1)),3), "\n chisq_val = ",round(1-Pearson_chi,3)))
カイ二乗値がPearson_chiに入っている
p値が0.95より小さければ棄却されない
棄却されない、とは
「有意水準5%（p値=0.95）において、クリック数がwebサイトに対して関係していない(結果からは関係していると言えない)」
ということ。
1-pchisq(Pearson_chi,1)
0.736
棄却されなかった
結果が正しかったのか組み込みの検定もしてみる
chisq.test(df,correct = FALSE)
correctをFにしておかないとYatesの補正が入ってしまうので注意
Pearson's Chi-squared test
data: df
X-squared = 0.11292, df = 1, p-value = 0.7368
計算の結果と検定の結果が一致していることを確認できた。
本当にカイ二乗分布に従うのかシミュレーションしてみる
とにかく二項分布発生させて複数回ぶん回す
chi_box = NULL
Pearson_chi = NULL
for(i in 1:10000){
click = rbinom(2,1000,0.8)
not_click = 1000-click
df = data.frame(click=click,not_click=not_click)
expected_frequency = colSums(df) * ( rowSums(df) / sum(rowSums(df)) )
col_1 = (df[1] - expected_frequency[1])^2 / expected_frequency[1]
col_2 = (df[2] - expected_frequency[2])^2 / expected_frequency[2]
Pearson_chi <- sum(col_1[1,1] + col_1[2,1] + col_2[1,1] + col_2[2,1] )
chi_box = c(chi_box, Pearson_chi)
}
h<-hist(chi_box, breaks=100)
dev.off()
plot(h$breaks[1:length(h$breaks)-1],h$density,type="n")
lines(h$breaks[1:length(h$breaks)-1],h$density)
d_chi = dchisq(h$breaks[1:length(h$breaks)-1],1)
lines(h$breaks[1:length(h$breaks)-1],d_chi, col="red")
histの密度と理論的な確率密度の値が一致するか描画
一致してるっぽいぞー
ピアソンではない方のカイ二乗値にできないか？
まずピアソンのカイ二乗値の計算でなく、
\frac{(X - np)}{\sqrt{np(1-p)}} = N(0,1) ??? \\
で計算した値は
「標準正規分布になるのか」
その後「標準正規分布の二乗がカイ二条分布に従うか」
を確認する
n=1000
norm_box = NULL
chi_box = NULL
for(i in 1:10000){
click = rbinom(2,1000,0.8)
not_click = 1000-click
df = data.frame(click=click,not_click=not_click)
expected_frequency = colSums(df) * ( rowSums(df) / sum(rowSums(df)) )
prob_binom_norm <- expected_frequency / sum(expected_frequency)
p = prob_binom_norm[1]
#極限としての正規分布近似の期待値と分散
mean_binom = n*p
var_binom =n*p*(1-p)
col_1 = (df[1] - mean_binom)^2 / var_binom
#col_2 = (df[2] - mean_binom)^2 / var_binom
setu_norm <- (df[1] - mean_binom) / sqrt(var_binom)
norm_box = c(norm_box, setu_norm[1,1],setu_norm[2,1])
setu_chi <- c(col_1[1,1] + col_1[2,1])
chi_box = c(chi_box, setu_chi)
}
par(mfrow=c(1,2))
hist(norm_box,breaks = 200,xlim=c(-5,5))
curve(dnorm(x,0,1),xlim=c(-5,5))
漸近的に標準正規分布に従っているといえそう・・・？
でも間違ってはいけないので詳しく見てみる
breaks=30
breaks=50
breaks=100
100にすると中央がくぼんでいるようにみえる
これは「近似」の限界なのか、何かミスってるのか？
自由度は1のはず
chisq_box <- chi_box^2
h<-hist(chisq_box, breaks=1000)
dev.off()
plot(h$breaks[1:length(h$breaks)-1],h$density,type="n",xlim=c(0,10))
lines(h$breaks[1:length(h$breaks)-1],h$density)
d_chi = dchisq(seq(0,20,0.1),1)
lines(seq(0,20,0.1),d_chi, col="blue")
これは似ているといっていいものか・・・？
まあ似ているしヨシっ!!
(ヨシじゃなかったら教えて)
この状態ならz検定や普通のカイ二乗検定もつかえるのでは？？
(使っちゃまずかったら教えて)
シミュレーション計算が合っているのかは別としても、理論側を調べた限り標準正規分布には近似できるらしい
ド・モアブル–ラプラスの定理
標本サイズn,確率pとする二項分布B(n,p)に従う確率変数Xは、
中心極限定理より
nが十分に大きい時、観測値Xは正規分布N(np,np(1-p))に近似する。
言い換えると
nが十分に大きい時、標本の比率X/nの従う分布がN(p,p(1-p)/n)の正規分布に従う。
(nをかけたら期待値np,分散np(1-p)になる)
標本の値を標準化することで
\frac{X-np}{\sqrt{np(1-p)}}
は標準正規分布となる。
なぜ「ピアソンのカイ二乗値」は分散でなく期待度数を分母に使ってカイ二乗分布に従うのか？を考える
標準正規分布になりそうなのは中心極限定理からわかった。
この問題の解法を書いている記事が全然見つからない
唯一阪大の先生の資料があったが線形代数が出てくるわ、脱字か故意かわからない部分が見られたりして理解はできなかった。
ので、今回の2×2のクロス集計表の「ピアソンのカイ二乗値」を変形・計算して、標準正規分布の二乗に置き換えることができれば万歳。
計算簡単のため、
・パラメータpはwebサイトA,Bのどちらでも同じである
・故にクリック数も同じであった
としておく。
click not_click \ 合計
A 800 200 \ 1000
B 800 200 \ 1000
-----------------------\
合計 1600 400
すべて記号に置き換える
全く同じだと面倒なので少しばかりズラす
click not_click 合計
A X+ε n-X-ε n
B X-ε n-X+ε n
ここで「ピアソンのカイ二乗値」は
\frac{(各セルの中の値 - 期待度数np)^2}{期待度数np}
をすべてのセルで求めて、すべてを合計するのだった。
\frac{X+ε - np}{np} + \frac{n-X-ε - n(1-p)}{n(1-p)} +\frac{X-ε - np}{np} + \frac{n-X+ε - n(1-p)}{n(1-p)}
簡単にして計算していく
\frac{(X - np)^2}{np} + \frac{(n-X - n(1-p))^2}{n(1-p)}
\frac{(1-p)*(X - np)^2}{np(1-p)} + \frac{p*(n-X - n(1-p))^2}{np(1-p)} \\
= \frac{(1-p)*(X - np)^2}{np(1-p)} + \frac{p*(X-np)^2}{np(1-p)} \\
= \frac{(X - np)^2}{np(1-p)} \\
4つのクロス集計表の値から1つの標準正規分布の二乗の形が得られた
期待値や分散としての役割を正規分布と置き替えても一致しそう。
\frac{(X - np)^2}{np(1-p)} = \frac{(X-\mu)^2}{\sigma^2} \\
故に自由度1のカイ二条分布に従いそうだな、って思うことができました(小並感)
これであってんの？
以上!!
参考
英語版
https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem
日本語版
https://k-san.link/binomial-normal/