library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

2.確率の復習

2.2 条件付き期待値

相関関係のある乱数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)
## [1] 0.713894
plot(y1, y2)

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

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

以下の式Aが成り立つことをシミュレーションにより確認する \[ Eq.A: \ E[y_2|y_1=0]= \int_{i=1}^{N}y_{2i}*p(y_{2i}|y_1=0)dy_2=\sum_{i=1}^{N}y_{2i}*Pr(y_{2i}|y_1=0) \]

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

#E[y2|y1=0]
df %>% 
  filter(y1==0) %>% 
  summarise(mean(y2))
##    mean(y2)
## 1 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
## # A tibble: 1 × 1
##   `sum(p * y2)`
##           <dbl>
## 1         0.307

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

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

以下の式Bが成り立つことをシミュレーションにより確認する \[ Eq.B: \ E[y_2]= \int_{i=1}^{N}E[y_2|y_{1i}]p(y_{1i})dy_1=\sum_{i=1}^{N}E[y_2|y_{1i}]*Pr(y_{1i}) \]

Eq.Bの左辺：E[y2]を求めてみる

#E[y2]
df %>% 
  summarise(mean(y2))
##   mean(y2)
## 1   -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)
## # A tibble: 1 × 1
##   `sum(e2 * p1)`
##            <dbl>
## 1       -0.00600

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

2.3.4 カイ二乗分布

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

#標準正規分布
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)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#実装済みのカイ二乗分布関数を使用
x <- c(1:50)
ggplot()+
  geom_line(aes(x = x,y=dchisq(x,1)))+
  geom_line(aes(x=x,y=dchisq(x,2)),colour="red")+
  geom_line(aes(x=x,y=dchisq(x,3)),colour="blue")+
  geom_line(aes(x=x,y=dchisq(x,4)),colour="green")+
  geom_line(aes(x=x,y=dchisq(x,20)),colour="purple")+
  ylim(0,0.3)
## Warning: Removed 1 row containing missing values (`geom_line()`).

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

#任意の自由度のカイ二乗分布のデータを作成する関数
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()

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

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

F分布を導出する

\[ F=((χ_1^2)⁄p)/((χ_2^2)⁄q) \]

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

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()