axjack's blog

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

標準誤差は統計量の標準偏差であることを確かめる

やること

離散型確率分布

X p
17 0.3
-1 0.7

に於いて、大きさn=100の標本を抽出する時、

  • 和S の標本分布について、Sの平均、標準誤差
  • 平均M の標本分布について、Mの平均、標準誤差

をそれぞれ

によりそれぞれ求める。

計算:離散型確率分布から分かること

Xの期待値・分散・標準偏差はそれぞれ、

  •  E(X) = \sum X_i p_i = 17\times0.3 + (-1) \times 0.7 = 4.4
  •  V(X)  = {SD(X)}^{2} \approx {8.249}^{2} \approx 68.04
  • SD(X) = |X_1 - X_2|\sqrt{p(1-p)} = |17 - (-1) |\sqrt{0.3\times0.7} \approx 8.249

であるから、サンプルサイズn=100における、

  • 和S =\sum{X_i} の標本分布について

    • X_i X_i \sim X として、
    • Sの期待値:=  E(S) = E(\sum X_i ) = E(nX_i) = nE(X) = 100 \times 4.4 = 440
    • Sの標準誤差:=  SE(S) = \sqrt{ V(S) } = \sqrt{ V(\sum{X_i} ) }  = \sqrt{ \sum{V(X_i)}  } = \sqrt{  \sum{V(X)} }  = \sqrt{ nV(X)} = \sqrt{n}\sqrt{V(X)}  = \sqrt{n} SD(X) = 10\times8.249 = 82.49
  • 平均M =\frac{1}{n}\sum{X_i} の標本分布について

    • X_i X_i \sim X として、
    • Mの期待値:=  E(M) = E(\frac{1}{n}\sum X_i ) = E(X) = 4.4
    • Mの標準誤差:= SE(M) = \sqrt{ V(M) } = \sqrt{  V(  \frac{1}{n}  \sum{X_i} ) } = \sqrt{  \frac{1}{n^{2} }  \sum{ V(X_i) }   }  = \frac{1}{n}  \sqrt{  \sum{ V(X)}  }  = \frac{1}{n} \sqrt{n}\sqrt{V(X)}  = \frac{1}{\sqrt{n}} SD(X) = \frac{8.249}{\sqrt{100}} = 0.8249

である。

モンテカルロ シミュレーション

共通の設定

X <- c(17,-1) #確率変数
p <- c(0.3,0.7) #確率
n <- 100 #サンプルサイズ
B <- 10000 #サンプル"数"

和の分布

# 和の分布 ####
S <- replicate(B,{
  sum(sample(X,size=n,replace=TRUE,prob=p))
})

# 結果
c(平均値 = mean(S), 標準誤差 = sd(S))
> c(平均値 = mean(S), 標準誤差 = sd(S))
   平均値  標準誤差 
438.75620  82.56416 

平均値の分布

# 平均値の分布
M <- replicate(B,{
  mean(sample(X,size=n,replace=TRUE,prob=p)) 
})

# 結果
c(平均値 = mean(M), 標準誤差 = sd(M) )
> c(平均値 = mean(M), 標準誤差 = sd(M) )
   平均値  標準誤差 
4.4000180 0.8259942 

おわりに

やってることは例によってChapter 13 Probability | Introduction to Data Scienceで学んだことです。標準誤差は統計量の標準偏差であるってのは数式だけだといまいちピンとこないのですが、このようにRでシミュレーションするとだいたい同じになるので腑に落ちました。replicate関数の第一引数にサンプル"数"を、第二引数には統計量(の式)を記述することによって、統計量の分布をシミュレーションすることができます。便利です。

おまけ: replicate関数の練習

(これを最初に書くべきだったのでは?)サイコロを2回投げて出た目の和が7となる確率、をRで求める。

共通

saikoro <- 1:6 # サイコロ
B <- 10000 #サンプル"数"

組み合わせで割り算する素朴な求め方

# サイコロを2回投げた時の組み合わせ
saikorokumi <- expand.grid(saikoro,saikoro)

# 和が7になる組み合わせ ÷ サイコロを2回投げた時の組み合わせの総数
sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi)
> sum( apply(saikorokumi,1,sum) == 7 ) / nrow(saikorokumi)
[1] 0.1666667

2回投げて出た目の和の分布

# サイコロを2回投げて出た目の和の分布
kiroku <- replicate(B,{
  sum(sample(saikoro,2,replace=TRUE))
})
# 和が7となる確率
mean(kiroku == 7)
> mean(kiroku == 7)
[1] 0.1659

2回投げて出た目の和が7となるか否かの分布

# サイコロを2回投げて出た目の和が7となるか否かの分布
kiroku2 <- replicate(B,{
  sum(sample(saikoro,2,replace=TRUE)) == 7 
})
# 和が7となる確率
mean(kiroku2)
> mean(kiroku2)
[1] 0.1684
axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.