2008年06月15日

ブートストラップとジャックナイフ

うーん。分かってきたような分からないような。

以前ブートストラップに関してちょっとだけ書いた。 あのときから別に大して詳しくなったわけではないが、講義の関係でクラスタ解析なんてものを少し勉強している。 DNAの変異だとかmRNAの発現の程度だとかいった情報から種や遺伝子のクラスタを作るときは、 その配列をブートストラップ法によって複製し、系統樹が支持される確率なんてのを計算するのが今では当たり前になっているそうだ。 この辺詳しく解説したらプレゼンの時間が稼げるかなとか企んでいる次第で。

まあそれはそれとして、ちょっとジャックナイフ法とブートストラップ法について勉強したのでメモ。 なんか前と同じような内容になりそうだ。まあいいや。

1.推定値

母集団の分布というのは、平均だとか分散だとかいった特性値によって特徴付けられる。もちろん普通はそれを直接に知ることはできない。 母集団の全てのデータを得ることは大抵の場合不可能だろう。そこで、 多くの場合は母集団からサンプリングして得られた標本より母集団の特性値を推定する。

母集団の平均値が知りたければ標本の平均値を計算するし、母集団の分散が知りたければ標本の分散を計算するだろう(もちろん、 分散の場合は標本分散ではなく普遍分散が推定量としては適切である)。

ここで問題となるのは、これらの推定値をどの程度信頼してもよいのか、つまり、これらの推定値はどれだけばらつきうるのか、 ということだろう。

平均値については前回ネタにしたので、今回は分散を扱うことにしよう。

2.カイ二乗分布を用いた分散の区間推定

詳細は面倒なので省くが、分散の区間推定はカイ二乗分布を用いて行われる。カイ二乗分布を利用した場合の95%信頼区間は、 データの平方和をS、自由度をfとして

20080615eq1

で表される。分母はそれぞれ自由度fのカイ二乗分布における2.5%点、97.5%点の値である。

しかし、この場合は母集団の正規分布を仮定している。母集団が正規分布しているとはいえない場合、 導かれる区間は適切なものとならない場合がある。

3.ブートストラップ法を用いた分散の区間推定

ブートストラップ法を用いる場合は母集団に正規分布のような分布を仮定しない。分散の区間は、 サンプリングからの複製によってその場で作り出される。詳細は前回解説したので省略。

4.ジャックナイフ法による分散の誤差推定

ジャックナイフ法とは1940年代にM. Quenouilleにより考案され、1950年代にJhon Tukeyにより命名された手法である。ブートストラップ法と異なるのは、 リサンプリングにおいていくつかのデータを抜いた状態のサンプルを生成するということ。狭義には1つのデータが除かれ、 重複サンプリングを行わない。そのため、ブートストラップ法よりも計算量は少ない。

4.1 ジャックナイフ法によるバイアス推定

推定値にはバイアスがかかる場合がある。たとえば、母分散の推定値として標本分散を採用することは妥当なことのように思えるが、 よく知られているように母分散の不偏推定量は不偏分散である。標本分散を母分散の推定値としてみた場合、 サンプル数に依存したバイアスがかかっているということになる。

n個のサンプル全てを用いて計算された推定量をθとする。n個のサンプルから、 i番目のデータを除いて作られた複製より計算された推定量をθiとする。その平均をθmeanとする。このとき、 ジャックナイフ法によりバイアスを修正された推定量は次式で計算される。

20080615eq2

たとえばこの式へ標本分散を代入して計算してみると、不偏分散となる。

4.2 ジャックナイフ法による誤差の推定

i番目のデータを除いて計算された推定値をθiとして、θの分散は

20080615eq3

によって計算される。

分散が分かるんだから信頼区間も分かるんだろうけどちょっと求め方がよくわからない。 θは中心極限定理から正規分布すると考えて良いんだろうか。

5. Rによる計算

ともかく、Rを用いて確認してみよう。

まず母集団の用意。

x <- 1:10

1〜10までの値が等しい確率で出現する離散分布と考えていい。分散は8.25、平均は5.5の一様分布となる。

ここから20個のデータからなる標本を抽出する。

#サンプルの生成
x.samp <- sample(x, 20, replace=T)

そして、カイ二乗分布を利用した場合とブートストラップ法による場合の区間推定。

#カイ二乗分布を利用した区間推定
19*var(x.samp)/qchisq(c(0.975,0.025),19)
#bootstrap法による区間推定
x.var.samp <- numeric(2000)
for(i in 1:2000){
 x.b <- sample(x.samp,20,replace=T)
 x.var.samp[i] <- var(x.b)
 }
quantile(x.var.samp,c(0.025,0.975))

次に、ジャックナイフ法によるバイアス修正と、「分散の分散」の推定。

#ジャックナイフ・バイアス修正済み推定量
x.j.var <- numeric(20)
for(i in 1:20){
 x.j.samp <- x.samp[-i]
 x.j.var[i] <- sum((x.j.samp-mean(x.j.samp))^2)/19
 }
sum((x.samp-mean(x.samp))^2)-19*mean(x.j.var)
#ジャックナイフ・分散推定
(19/20)*sum((x.j.var-mean(x.j.var))^2)

手元で一回実行してみたところ、

  • サンプルの不偏分散: 7.671053
  • ジャックナイフ・バイアス修正済み標本分散: 7.671053
  • カイ二乗分布による信頼区間(95%): 4.43652 - 16.36442
  • ブートストラップ法による信頼区間(95%): 4.365592 - 10.450526
  • ジャックナイフ法による分散の分散の推定: 2.550627
  • ブートストラップ法による分散の分散: 2.371628

となった。まず、サンプルの不偏分散とジャックナイフ・バイアス修正済み標本分散が同じになるのは当然。計算の中身が同じなので。

カイ二乗分布による区間はブートストラップによるものより広く出ている。カイ二乗分布による区間はほぼ100% の確率で母分散である8.25を含んだ。母集団が正規分布でない影響だろうか。一方、 ブートストラップ法による区間は9割前後の確率で母分散を含んだ。2000回程度の試行をしてみたが、 母分散が信頼区間に含まれる確率がどうも95%よりも低い気がする。なんでだろ。

ジャックナイフ法による分散の分散の推定とブートストラップ法による分散の分散は大体同じになった。うーん。 両者がそれほど離れることはないし、これだとジャックナイフ法で十分じゃないかということになるなぁ。分散だから? ブートストラップは適用範囲が広いとかいう話も見たなー。

まーとにかく今日は眠いので寝よう。参考文献はまたあとで追加。追加しました。

参考

posted by Rion778 at 03:40 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2008年06月09日

円周率と根長

格子法 (根長を測定するほうの)ってのは要するにBuffonの針と同じ問題、というかむしろBuffonの針の問題を応用しているんだろうか。 論文はよく確認していなかったけど、緒言とか参考文献にそれらしいものが書いてあったのかもしれない。

そんなことを思いついたのだが、もう眠いので寝る。時間があったらいつか調べる。

posted by Rion778 at 23:53 | Comment(0) | TrackBack(0) | diary | このブログの読者になる | 更新情報をチェックする

円分方程式

20080609eq1

以前、Rを使って正n角形を描く関数を定義したとき、 上記の円分方程式というものを利用した。この方程式の根は、複素平面上の単位円の円周をn等分する、というものであった。

この方程式の根は、つまり1のn乗根ということで、代数学の基本定理よりそれはn個存在する。

たとえば、次の3次方程式の場合ならば、

20080609eq2

この式の根は次の3つである。

20080609eq3

そしてこれを複素平面上にプロットし、それぞれ直線で繋げば正三角形が描かれる。

tri

そして、この根はnがいくつだろうと次の式によってすべて求められる。

20080609eq10

以下、ちっとだけ詳しく解説。

まず、この複素平面というやつだが、「a+bi」という形で表される複素数の、aを横軸の「実軸」に、bを縦軸の「虚軸」 に対応させることで平面上に複素数を表現している。平面上のひとつの点が、ひとつの複素数に対応する。

そして、この平面上の単位円ということだから、単位円周上の複素数のaとbはすべて

20080609eq4

を満足するということになる。左辺の式で表現されるものは複素数の絶対値と呼ばれる。つまり、 絶対値が1の複素数の集合が複素平面の単位円を形成する。また、複素数の絶対値とは極座標表現における「長さ」 であることもわかる。

ここで、原点から複素数へ引かれた直線と実軸とがなす角、複素数の絶対値の2つの情報があれば、三角関数を用いてaとbを表現できる。 絶対値をr、なす角をθとし、aとbはそれぞれ次のように表現できる。

20080609eq5

ここで角θは偏角と呼ばれる。ある複素数zの偏角である、ということを意識したい場合は、「θ= arg z」などと書く。複素数zは、aとbの代わりに偏角と絶対値を用いて次のように表現できる。

20080609eq6

このような複素数の表現方法を極形式と呼ぶ。極形式を用いて、複素数同士の積を考察する。

20080609eq7

最後の変形には三角関数の加法定理を用いている。この変形から明らかなように、 複素数同士の積の絶対値はもとの複素数の絶対値の積であり、 偏角はもとの複素数の偏角の和となる。

単位円の場合、絶対値は1であるため、偏角θを持つ単位円上の複素数を乗ずるということは、 複素平面上におけるθ分の回転を意味する

話を1のn乗根へ戻そう。1のn乗根ということは、n回乗ずると1になるような数、ということである。 先の複素数の積に関する考察を利用して言い換えれば、 n回回転させると角度が0(=2π)となるような角θを偏角に持つ絶対値1の複素数、ということである。 単位円上の複素数をn回乗ずると、偏角はn倍される

さて、n倍すると2πになるような角はθ=0を含め、n個存在する。それらはすべて2π/nの整数倍である(よく確認してほしい)。 このことに注意すれば、

20080609eq1

の根は

20080609eq8

で表されるn個の複素数であるということがわかる。

ここに、オイラーの公式を適用して…ってブログ内に記事がないじゃないか。とっくに書いてるものだと。 仕方ないので今日のところは天下り式に次の式を。(オイラーの公式について簡潔にまとまっていて分かり易い解説:オイラーの公式[物理のかぎしっぽ])

20080609eq9

これがオイラーの公式と呼ばれるもので、θにπを代入したりすると大変なことになったりするが今は触れない。 オイラーの公式と上の式を見比べると、1のn乗根は次の式で求められることがわかる。

20080609eq10

まあ書き方がシンプルになっただけで中身は変わっていない。オイラーの公式を使うと、 極形式の複素数を指数の形でひとつにまとめられる、ということである。

それで、以前Rを使って正n角形を書く関数を定義したときはこの最後の式だけを使ったわけだが、 当時は今ひとつこの式の意味がわからないまま使っていた気がする。

さっき寝ようと思って布団の中で本を読んでいたら1のべき根の話が書いてあって、 1のn乗根を求める上の式の意味がなんとなく分かってきたので、いろいろなことを棚上げして夜更かしして記事にまとめた次第。

参考文献

posted by Rion778 at 02:43 | Comment(0) | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。