標本分布の平均と分散について、実際にプログラムを使って検証してみます。
標本分布の平均と分散
X を、ある分布(distribution)に従う確率変数(random variable)とします。 X1,⋯,Xm を独立な X の標本とし、 Y を次のように定義します。
Y=1mm∑k=1Xk
このとき、 E(Y), V(Y) は次のようになります。
E(Y)=E(X)V(Y)=frac1mE(X)
説明
標本分布の平均と分散の式の算出をみておきましょう。
Xi は独立なので E(Y)=1m∑mk=1E(Xk)=1mE(X) です。
分散についても 2つの確率変数 A, B が独立ならば V(A+B)=V(A)+V(B) となる。 これより V(Y)=mV(1mX)=1mV(X) 。
2項分布を使った実験
本当にそうなるのか、2項分布を例に、コンピュータを使って実験してみます。
X∼textrmBinomial(n,p) とすれば E(X)=np, V(X)=np(1−p) です。
Y を m個の独立な標本の平均値とします。 このとき E(Y)=1mnp),\(V(Y)=1mnp(1−p) です。
ではここで n=3, p=0.5, m=10 として計算してみます。 E(X)=E(Y)=0.5, V(X)=0.75, V(Y)=0.075 となります。
Kotlin によるコード
まず n
, p
, m
をそれぞれ定義しておきます。
|
val n = 3 val p = 0.5 val m = 10 |
それぞれのベルヌーイ試行から、 X の値をシミュレーションするメソッドを定義します。 Math.random
を使って 0以上1未満の乱数を生成し、 p
と比較して、 p以上なら 1を足す操作を n 回 行なっています。
|
fun exec(n: Int, p: Double): Int { var result = 0 for (i in 1..n) if (Math.random() >= p) result++; return result } |
この試行の結果を m回取得して平均を計算するメソッドを定義します。 これは Y の値の計算にあたります。
|
fun exec(m: Int, n: Int, p: Double): Double { var result = 0.0 for (i in 1..m) result += exec(n, p) return result / m } |
(Y) を10000回計算して、平均と分散を計算します。
|
val timeCount = 100000 var mean = 0.0 var variance = 0.0 repeat(timeCount) { val result = exec(m, n, p) mean += result variance += Math.pow(exec(m, n, p) - 1.5, 2.0) } // np = 1.5 println(mean / timeCount) // np(1-p)/m = 0.075 println(variance / timeCount) |
以上をまとめると次のようになります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
|
package com.improve_future.sample object Main { @JvmStatic fun main(args: Array<String>) { val n = 3 val p = 0.5 val m = 10 val timeCount = 100000 var mean = 0.0 var variance = 0.0 repeat(timeCount) { val result = exec(m, n, p) mean += result variance += Math.pow(exec(m, n, p) - 1.5, 2.0) } // np = 1.5 print("Expectation: ") println(mean / timeCount) // np(1-p)/m = 0.075 print("Variance : ") println(variance / timeCount) } fun exec(n: Int, p: Double): Int { var result = 0 for (i in 1..n) if (Math.random() >= p) result++; return result } fun exec(m: Int, n: Int, p: Double): Double { var result = 0.0 for (i in 1..m) result += exec(n, p) return result / m } } |
このコードの出力は次のようになり、 計算通りの E(Y), V(Y) の値になっていることが確認できます。
|
Expectation: 1.4999001999999781 Variance : 0.00756705900000417 |
A Life Summary of an Gypsy