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

Table of Contents

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

標本分布の平均と分散

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

Y = \frac{1}{m} \sum_{k = 1}^{m} X_{k}

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

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

\begin{array}{rcl} E (Y) & = & E (X) \\ V (Y) & = & f r a c 1 m E (X) \end{array}

$\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>) {