axjack's blog

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

不偏分散が不偏であることを確認する

元ネタ

統計の基礎母集団からランダムに取り出した標本でなくても,例えば10個の固定した値を4個と6個に分ける場合,両者を比較するには何で割った分散を使うかという問題についても,一貫して n − 1で割るほうが比較としては正しくなります: を別データに差し替えて確認する。

何を確認するのか

大きさN=9 の xdataを母集団とみなし、大きさn = {2, ... ,9} の標本を抽出する時の、それぞれの不偏分散を確認する。

コード

#データを用意し
xdata <- c(11,22,33,44,55,66,77,88,99)

#不偏分散を自力で計算
huhenbunsan <- (1/( length(xdata) - 1)) * sum( ( xdata - mean(xdata) )^2)

#あらかじめ用意されている不偏分散の関数
var(xdata)
#> var(xdata)
#[1] 907.5

#答え合わせ
var(xdata) - huhenbunsan
#> var(xdata) - huhenbunsan
#[1] 0

#xdataから任意の2つの組み合わせで不偏分散を取り
#それらの平均値を求める
mean( combn(xdata,2,var) )
#> mean( combn(xdata,2,var) )
#[1] 907.5

#xdataから任意の3つの略
mean( combn(xdata,3,var) )
#> mean( combn(xdata,3,var) )
#[1] 907.5

#xdataから任意の4つの略
mean( combn(xdata,4,var) )
#> mean( combn(xdata,4,var) )
#[1] 907.5

#xdataから任意の5つの略
mean( combn(xdata,5,var) )
#> mean( combn(xdata,5,var) )
#[1] 907.5

#xdataから任意の6つの略
mean( combn(xdata,6,var) )
#> mean( combn(xdata,6,var) )
#[1] 907.5

#xdataから任意の7つの略
mean( combn(xdata,7,var) )
#> mean( combn(xdata,7,var) )
#[1] 907.5

#xdataから任意の8つの略
mean( combn(xdata,8,var) )
#> mean( combn(xdata,8,var) )
#[1] 907.5

#xdataから任意の9つの略
mean( combn(xdata,9,var) )
#> mean( combn(xdata,9,var) )
#[1] 907.5

おわりに

感想

標本の不偏分散の平均値が、母集団の不偏分散と全部一致していて気持ちいい(?)です。

combnについて

combn {utils}    R Documentation
Generate All Combinations of n Elements, Taken m at a Time

Description

Generate all combinations of the elements of x taken m at a time. 
If x is a positive integer, returns all combinations of the elements of 
seq(x) taken m at a time. If argument FUN is not NULL, 
applies a function given by the argument to each point. If simplify is FALSE, 
returns a list; otherwise returns an array, typically a matrix. ... 
are passed unchanged to the FUN function, if specified.