【統計】分散について【R】

Pocket

こんにちは、まんぼうです。

今回も統計に関する記事です。

タイトルにある通り分散についてまとめていきます(この体の記事は私の備忘録も兼ねてます。)

はじめに

分散とは一言で言うとデータの散らばり具合を示すもの。変動を評価するもの。

平方和

平方和は、式で表すと

こう。

Tex記法で書こうと思ったのですが、上手くいかなかったので他のところで書いたもののスクリーンショットです。

ちなみにTex記法で書くと

\sum_{i=1}^{n}(yi-\hat{y})^2

こう。

2乗しているのはしないとどんなデータでも0になるからです。

絶対値を取るという方法も考えられますが、それは計算が難しくなるみたい。

なので残差([$ (yi-\hat{y})])のことを二乗してからその総和を取るのが良いらしいです。

自由度

一般的に自由度の正式な定義は、「自由度は n-p である。ただしnは標本数であり、pは推定された母数の数である。」

例えば標本数が5つのデータがあるとしてその平均が求まっているとした場合、n=5でp=1であるため、自由度は４であると言える。

分散

ここまでは前置きで、ここから分散について。

分散は s^2 で表され、分散＝平方和/自由度、で求めることが出来ます。

つまり式で表すと

Tex記法で書くと

[$ s^2=\frac{\sum_{i=1}^{n}(yi-\hat{y})^2}{n-p}]

一般的にp=1なので

と表されるし、高校の時の数学ではこう習った覚えがある。

これもまたTex記法で書くと

[$ s^2=\frac{\sum_{i=1}^{n}(yi-\hat{y})^2}{n-1}]

厳密には分母の-1は自由度を考えている。

これをRで書くと

sum((x-mean(x))^2)/(length(x)-1)

となる。

この組み込み関数は

var(data)

分散が同じで平均が異なる場合

2つの母平均を比較する標準的な検定(スチューデントのt検定など)や3つ以上の母平均を比較する分散分析が適切

平均が同じで分散が異なる場合

この場合、というより分散が異なる場合平均を比較して推測すべきではない。

その分散の違いは有意であるかどうかを調べるためにはF検定を行う。まずは分散が小さいの方の値で大きい方の値を割る。

var <- var(data1)/var(data2)

(var(data1)>var(data2))とする

この時に2つの母分散が真に等しいという仮説の元で、F比が偶然によってこの値よりも大きくなる確率を調べる。

これにF分布の累積確率分布が必要。これはRではpf関数と呼ばれ、分散比(上のvar)、分子の自由度(f1とする)、分母の自由度(f2とする)を必要とする。

2*(1-pf(var, f1, f2))

２掛けてあるのは両側検定というものを実行するため(data1とdata2でどちらの方が大きい分散を持っているのかが分かっていないため)。

ちなみにこれは、Rの組み込み関数でも出来る。

var.test(data1, data2)

分散と標本数

plot(c(0,32),c(0,15),type="n",xlab="Sample size", ylab="Variance")

for (df in seq(3,31,2)) {
    for (i in 1:30) {
        x<-rnorm(df,mean=10,sd=2)
        points(df, var(x))
    }
}

このコードの実行結果

標本数が減少すると推定値の範囲が増加している。標準偏差(sd)はsd=2とされているため分散は4で一定である。

ただグラフを見れば明らかに標本数が大きいほど推定値の範囲が減少し、標本数が減少すると推定値の範囲は増加する。

ちなみにこれは、サンプル数30が多いか少ないかの境目であることの根拠になっている。

分散の使用目的

・非信頼度を評価するため（例えば信頼区間の設定）

・仮説検定を行うため（例えばスチューデントのｔ検定）

非信頼度の指標

データの分散が増加するときに非信頼度は増加すべきであるから、非信頼度∝[$ s^2]となる。

一方標本数は増加すると非信頼度は減少すべきであるため、非信頼度∝ s^2/n となる。

ただしこの単位は何かの2乗になる。これは指標の対象となる次元とは同じであることが望ましいため好ましくない。そのため平方根を取る。

Tex記法で書くと[$ s.e(\bar{y})=\sqrt{\frac{s^2}{n}}]

こうなる。これは”標準誤差”という。ちなみにこの前の平方根を取っていないものは”平均の標準誤差”という。

信頼区間

信頼区間とは標本に基づいて作られる区間で、それらが母平均を含む可能性の高い区間のこと。

非信頼度が大きくなると、信頼区間が広くなる。つまり式に表すと

(信頼区間の幅)∝(非信頼度の指標)∝[$ \sqrt{\frac{s^2}{n}}]

となる。(分散が大きいと幅も大きくなるよねっていうこと)

確率点を求めるため、各水準の信頼区間に関連したスチューデントのｔ値はqt関数で求められる。これはｔ分布の分位点を計算する。

例えば95％の信頼区間を求めたいときはα=0.025(つまり、0.01*(100%-95%)/2←/2しているのは左右両方に信頼区間でない区間があるから)に関連したｔ値を利用する。左側は0.025の分位点、右側は0.975の分位点が必要になる。

qt(0.025,9)

qt(0.975,9)

この第1引数は確率で第2引数は自由度。

この第1引数の確率が大きくなると信頼度の幅も広くなるためｔ値も大きくなる。

小さな標本数の時の母平均の信頼区間に対する公式は

(信頼区間)=(標本平均)±(ｔ値)×(標準誤差)

となる。

ブーストラップ

信頼区間を計算する全く別の方法

data<-read.table("c:/~~/", header=T)
attach(data)
names(data)

plot(c(0,30),c(0,60),type="n",xlab="Sample size",ylab="Confidence interbal")
for (k in seq(6,30,3)){
    a<-numeric(10000)
    for (i in 1:10000) {
        a[i]<-mean(sample(values,k,replace=T))
    }
    points(c(k,k),quantile(a,c(0.025,0.975)),type="b")
}

信頼区間の幅はn=20くらいまでは急速に狭くなっていくがそのあとはゆっくりになる。

n=30で10000回のシミュレーションにも基づくブーストラップCIは次のように求める

quantile(a,c(0.025,0.975))

まんぼう放浪記