Rで条件付き期待値の導出確認、χ二乗分布、F分布の確認

一般化線形モデルの勉強の一環として、条件付き期待値と確率分布の性質についてRの乱数シミュレーションで確認してみる

条件付き期待値の性質

相関関係のある乱数y1,y2を生成

set.seed(1)
rho <- 0.7
x <- rnorm(n=1000,mean=0,sd=3)
e1 <- rnorm(n=1000,mean=0,sd=3)
e2 <- rnorm(n=1000,mean=0,sd=3)
y1 <- sqrt(rho)*x+sqrt(1-rho)*e1
y2 <- sqrt(rho)*x+sqrt(1-rho)*e2
cor(y1,y2)
plot(y1, y2)

df<-as.data.frame(cbind(y1,y2)) %>% 
  round(0)

条件付き期待値：y1=0のときのy2の期待値E[y2|y1=0]を求める

以下の式Aが成り立つことを確認する

Eq.Aの左辺：E[y2|y1=0]を求めてみる

#E[y2|y1=0]
df %>% 
  filter(y1==0) %>% 
  summarise(mean(y2))

0.3070866

Eq.Aの右辺を求める

#E[y2|y1=0]の推定値 Pr(y2|y1=0)*y2の合計
df %>% 
  filter(y1==0) %>% 
  select(y2) %>% 
  group_by(y2) %>%
  summarise(n=n()) %>% 
  ungroup() %>%
  mutate(total=sum(n),
         p=n/total) %>% #Pr(y2|y1=0)
  summarise(sum(p*y2)) #sum of Pr(y2|y1=0)*y2

0.3070866

上記２つの導出法においてy2の条件付き期待値E(y2|y1=0)は一致した

周辺期待値E[y2]を求める

以下の式Bが成り立つことを確認する Eq.Bの左辺：E[y2]を求めてみる

#E[y2]
df %>% 
  summarise(mean(y2))

-0.006

Eq.Bの右辺：sum of E[y2|y1i]*Pr(y1i)を求めてみる

#Pr(y1)
Pr.y1<-df %>% 
  group_by(y1) %>%
  summarise(n=n()) %>% 
  ungroup() %>%
  mutate(total=sum(n),
         p1=n/total) %>% 
  select(y1,p1)


df %>% 
  group_by(y1) %>% 
  summarise(e2=mean(y2)) %>% #E[y2|y1]
  left_join(Pr.y1,by="y1") %>% 
  mutate(e2*p1) %>% #E[y2|y1]*Pr(y1)
  summarise(sum(e2*p1)) #sum of E[y2|y1]*Pr(y1)

-0.006
上記２つの導出法においてy2の周辺期待値の結果は一致した。

標準正規分布からカイ二乗分布を導出する

#標準正規分布
rnorm(1000) %>% 
  hist()

#df=1のカイ二乗分布
rnorm(1000)^2 %>% 
  hist()

#df=2のカイ二乗分布
df2<-rnorm(1000)^2+rnorm(1000)^2
df2 %>% 
  hist()

#以下は間違っている；独立な確率変数の二乗和である必要があるため。以下は独立でない単一の確率変数の整数倍になっている
a<-2*rnorm(1000)^2
a %>% 
  hist()


#任意の自由度のカイ二乗分布をヒストグラムで作成する関数
f.chi<-
  function(df=1){
  replicate(1000,#下記操作を1000回繰り返し、ヒストグラムを作成する
          rnorm(df)^2 %>% sum()) %>%   #標準正規分布からdf個の要素を取り出し、二乗和を求める                                           
          as.data.frame() %>% 
          ggplot(aes(x=.))+
          geom_histogram()
  }

f.chi(1);f.chi(2);f.chi(5);f.chi(100)

カイ二乗分布は自由度というパラメータのみで規定されている。また、自由度を増やすと正規分布に近づくこともわかる。

カイ二乗分布の再生性を確認する

#任意の自由度のカイ二乗分布のデータを作成する関数
chi<-function(df1=1,df2=1){
  replicate(1000,#下記操作を1000回繰り返す
          rnorm(df1)^2 %>% sum()) #標準正規分布からdf個の要素を取り出し、二乗和を求める
}

#任意の自由度の２つのカイ二乗分布のデータを作成する関数
chid<-function(df1=1,df2=1){
          a<-rnorm(df1)^2 %>% sum()
          b<-rnorm(df2)^2 %>% sum()
          c<-a+b
}
chi2<-function(df1,df2){
replicate(1000,chid(df1,df2))
}

#df=8のカイ二乗分布
chi(8) %>% 
  hist()

#df=3, df=5のカイ二乗分布
chi2(3,5) %>% 
  hist()

自由度8のカイ二乗分布（左）と自由度3,5のカイ二乗分布を足し合わせた確率分布は形状がかなり類似していることがわかる。再生性をデータから再現できた。

F分布を標準正規分布→カイ二乗分布から作成する

自由度p,qのカイ二乗分布があったとき、F統計量=(χ1/p)/(χ2/q)は自由度(p,q)のF分布に従う

f<-function(df1=1,df2=1){
          a<-rnorm(df1)^2 %>% sum()
          b<-rnorm(df2)^2 %>% sum()
          c<-(a/df1)/(b/df2) #F分布の定義
}

f.dis<-function(df1,df2){
replicate(1000,f(df1,df2))
}

#自作関数のF分布(df=10, df=20)
f.dis(10,20) %>% 
  hist()

#R標準のF分布(df=10,df=20)
rf(1000,10,20) %>% 
  hist()

自由度(10,20)のF分布を標準正規分布→カイ二乗分布から導出することができた

数式だけ見ていても、確率分布を二乗して足すとか割るって一体どういうことなのかよくわかっていなかったが、実際にRで乱数から確率分布を導出してみると納得した。数式では確率分布のまま扱っているが、実際には確率分布から標本抽出された（確率分布に従う）実現値を二乗して、それを1000個とかたくさん無作為抽出してヒストグラムにしたとき、その標本集団がこういう確率分布に従うんだよ、っていうことを意味しているということが直感的に理解できた。やはり理論の勉強はシミュレーションを併用するとかなり学習効率が高い。