axjack's blog

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

「骨波田の藤」の「骨波田」を何と読むか

「こっぱた」もしくは「こっぱだ」と読むものだと思っていたらどうやら「こつはた」と読むのが正しいらしので、インターネットから調べられる範囲で調べてみることにした。

「こつはた」と読むケース

「こっぱた」もしくは「こっぱだ」と読むケース

二変量正規分布の条件付き期待値・条件付き分散

以下のような二変量正規分布を考える。

 
\begin{bmatrix}
   X \\
   Y
\end{bmatrix}
\sim
N(
\begin{bmatrix}
   \mu_{X} \\
   \mu_{Y}
\end{bmatrix}
,
\begin{bmatrix}
   \sigma_{X} & \sigma_{XY} \\
   \sigma_{XY} & \sigma_{Y}
\end{bmatrix}

)

ここで、


 Z = X + AY, A = - \sigma_{XY} {\sigma_{Y}} ^{-1}

というZについて、 Cov(Z, Y) を計算すると、


cov(Z,Y) \\
= cov(X+AY,Y) \\
= cov(X,Y) + cov(AY,Y) \\
= \sigma_{XY} + Acov(Y,Y) \\
= \sigma_{XY}  - \sigma_{XY} {\sigma_{Y}} ^{-1} \sigma_{Y} \\
= \sigma_{XY}  - \sigma_{XY} \\
= 0

となる。 ゆえにX,Y,Zは正規分布に従う確率変数であることから、ZとYは共分散が0⇔ZとYは独立となる。

次に、 E[ X | Y ] を計算する。ここで X = Z - AY を用いると、


 E[ X | Y ]  \\
= E[  Z - AY | Y ]  \\
= E[  Z | Y ] - E[ AY | Y ]  \\
= E[  Z  ]   - A E[ Y | Y ] \\
= E[  X + AY ]  - AY   \\
= E[  X ]  + AE[Y] - AY \\
= \mu_{X} + A\mu_{Y} - AY \\
= \mu_{X} - A(  Y - \mu_{Y}  )  \\
= \mu_{X} + \sigma_{XY} {  \sigma_{Y}  }^{-1}  (  Y - \mu_{Y}  )

となる。*1

最後に、  V[ X | Y ] を計算する。


V[ X | Y]  \\
= V[ Z - AY | Y ] \\
= V[ Z | Y ] + V[AY | Y] - 2Cov(Z,AY) \\
= V[ Z | Y ] + V[AY | Y] - 2ACov(Z,Y) \\
= V[ Z | Y ]  + 0 + 0\\
= V[ Z ]

より、*2


V[ Z ]   \\
= V[ X + AY ] \\
= V[X ]  + V[ AY ] + 2Cov(X,AY)   \\
= V[X ]  + A^2V[ Y ] + 2ACov(X,Y)   \\
= \sigma_{X} + A^2\sigma_{Y} + 2A\sigma_{XY} \\
= \sigma_{X} + (- \sigma_{XY} {\sigma_{Y}} ^{-1})^2 \sigma_{Y} + 2(- \sigma_{XY} {\sigma_{Y}} ^{-1}) \sigma_{XY} \\
= \sigma_{X} + {\sigma_{XY}}^2 {\sigma_{Y}}^{-1} - 2{\sigma_{XY}}^2 {\sigma_{Y}}^{-1} \\
= \sigma_{X} - {\sigma_{XY}}^2 {\sigma_{Y}}^{-1}

となる。

参考

stats.stackexchange.com

*1:途中、性質1:XとYが独立ならばE[X|Y] = E[X] と性質2:E[X|X] = X を用いた。

*2:条件付き分散の性質がよくわからない・・・。リンク参照。

正規分布と適合度検定

適合度検定 :: 株式会社アイスタット|統計分析研究所より「適合度の検定(正規性)の結果」をRにて計算してみる。

# パラメータ####
# 平均
m1  <- 64.5
# 標準偏差
sd1 <- 13.41
# 度数の総和
n   <- 40

# 関数 ####
# 区間a<x<bにおける標準正規分布に従うXの確率、を返す関数
f <- function(l,h){ pnorm(h,m1,sd1) - pnorm(l,m1,sd1) }

# データ
# 観測値:observed ####
c(2,4,7,13,10,3,1) -> obs

# 期待度数expected ####
n * c(
pnorm(40, m1, sd1)
,f(40,50)
,f(50,60)
,f(60,70)
,f(70,80)
,f(80,90)
,(1 - pnorm(90, m1, sd1) ) 
) -> expected

# 自由度 = 階級の個数 - 1 - 2
df1 <- 7-1-2

# χ²統計量
chi2 <- sum( (obs - expected)^2/expected )
print( chi2 )
#> print( chi2 )
#[1] 1.382639

# p値
1 - pchisq( chi2, df1) 
#> 1 - pchisq( chi2, df1) 
#[1] 0.8472068
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.