標本分布の平均と分散を検証してみる

標本分布の平均と分散について、実際にプログラムを使って検証してみます。

標本分布の平均と分散

\(X\) を、ある分布(distribution)に従う確率変数(random variable)とします。 \(X_1 , \cdots , X_m\) を独立な \(X\) の標本とし、 \(Y\) を次のように定義します。

\[ Y = \frac{1}{m} \sum _{k=1}^{m} X_k\]

このとき、 \(E(Y)\), \(V(Y)\) は次のようになります。

\begin{eqnarray*} E(Y) & = & E(X) \\ V(Y) & = & frac{1}{m} E(X) \end{eqnarray*}

説明

標本分布の平均と分散の式の算出をみておきましょう。

\(X_i\) は独立なので \(E(Y) = \frac{1}{m} \sum _{k=1}^{m} E(X_k) = \frac{1}{m} E(X)\) です。

分散についても 2つの確率変数 \(A\), \(B\) が独立ならば \(V(A+B) = V(A)+V(B)\) となる。これより \(V(Y) = m V \left( \frac{1}{m} X \right) = \frac{1}{m} V(X) \) 。

2項分布を使った実験

本当にそうなるのか、2項分布を例に、コンピュータを使って実験してみます。

\(X \sim textrm{Binomial}(n, p)\) とすれば \(E(X) = np\), \(V(X) = np(1-p)\) です。

\(Y\) を \(m\)個の独立な標本の平均値とします。このとき \(E(Y) = \frac{1}{m} np), \(V(Y) = \frac{1}{m} np(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

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() &gt;= p) result++;
return result
}

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
}

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)

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)

以上をまとめると次のようになります。

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
}
}

package com.improve_future.sample

object Main {

@JvmStatic

fun main(args: Array<String>) {