フィッシャーの検定とモンテカルロ法③

以下の記事からの続き。

abcxyzonetwothree.hatenablog.com

統計を使って何らかの事象を分析する統計学のエンドユーザーにとって、モンテカルロ法は「コンピューターでランダムな数字を発生させて、理論から確率を求めるのではなく、経験から確率を求めよう」というざっくりとした説明で、基本的には事足りると思われる。

このブログでも、ランダムに数を発生させて、色んな実験を行ってきたが、それもある意味モンテカルロ法だと思う。

最も分かりやすのが、サイコロの目の出る確率の例である。

目の前に歪んだサイコロがあるとする。この歪んだサイコロは6つの目を持つ。しかし、歪んでいるため、例えば1の出る目の確率はわからない。どのようにその確率を求めれば良いのか。

丁寧な人は、その歪み具合をいろいろ計算して理論的に確率を求めるかもしれない。それによって、綺麗に計算できたら良いが、いろいろ大変そうである。そこで、何度もサイコロを振って、確率を経験的に求めようとする発想があり得る。これがモンテカルロ法の基本的な発想だと思われる。

もう少し具体的な例として、モンテカルロ法によって円周率を求めるというものがある。色んな記事がすでに解説されているが、改めてここでも解説する。

当たり前だが、円の面積は以下の公式によって求められる。

$円の面積=r^2×π$

なお、ここでrはその円の半径である。また、中心角が90°の扇型の面積は、以下のようになる。

$中心角が90°の扇型の面積=円の面積×\frac{90}{360}=r^2×π×\frac{1}{4}$ （１式）

ここで、以下のように考える。

まず、縦1cm、横1cmの正方形を考える。その中に、半径が1cmの中心角が90°の扇型が入っている。すなわち、以下のような図を考える（画像の比で、正方形には見えませんが、ご了承ください）。

f:id:abcxyzonetwothree:20181125115508p:plain

この中にランダムに点を打つ。例えば、100個の点を打った時、扇型の中に打たれる点（青色）と、外に打たれる点（赤色）が出てくる。

f:id:abcxyzonetwothree:20181125115602p:plain

この時、青色の点の数と全ての点の数（赤色+青色）の比は、扇型の面積と正方形の面積の比であると言える。また、この図において、正方形の面積は1×1で1なので、ここでの扇型の面積は以下のように表現できる。

$中心角が90°の扇型の面積=\frac{青色の点の数}{全ての点の数}×正方形の面積 = \frac{青色の点の数}{全ての点の数}$ （２式）

ここで、この扇型の半径は1cmなので、１式と２式から、円周率は以下のように表現できる。

$r^2×π×\frac{1}{4} = \frac{青色の点の数}{全ての点の数}⇄π=4×\frac{青色の点の数}{全ての点の数}$

この式を用いて、ランダムに打った点の数を用いて、円周率の推定値を求めることができる。

これがモンテカルロ法による円周率の求め方である。直感的にもわかると思うが、最初のランダムに打つ点の数が上がれば上がるほど、その精度は良くなる。

実際に、上記の100個の点で試した場合の円周率の推定値は3.32であった。

1000個の場合は、以下の図のようになり、その推定値は3.164であった。

f:id:abcxyzonetwothree:20181125115748p:plain

10000個の場合は、以下の図のようになり、その推定値は3.1408であった。

f:id:abcxyzonetwothree:20181125115758p:plain

なお、そのコードは以下の通り。

#円の関数
f<-function(x){
  sqrt(1-x^2)
}

#x軸の設定
x<-seq(-1,1,0.01)

#扇型の図示
curve(f(x),ylim=c(0,1),xlim=c(0,1),xlab = "x",ylab="y")

n<-10000 #乱数の発生回数（ここを、100、1000にする）
a<-runif(n, min=0, max=1)
b<-runif(n, min=0, max=1)

#点の図示
par(new=T)
plot(a,b,col=ifelse(b>sqrt(1-a^2),"red","blue"),ylim=c(0,1),xlim=c(0,1),,xlab = "",ylab="")

#円周率の計算
count<-sum(ifelse(b>sqrt(1-a^2),0,1))
(count*4)/n

それでは、フィッシャーの検定の時にどのように使われているのか。いろいろ、文献を漁ったが、非常に難しく正確なことは述べられないが、ざっくりとまとめると以下のような説明となる。

前回の記事でも説明した通り、フィッシャーの検定では、行列の事象が独立という条件のもと、周辺度数を固定した上で、あり得る全てのマトリクッスを求め、それぞれにおいてその出現確率を計算する。前回の記事では、恣意的に簡単なマトリックスを用意したが、セル度数が多くなったり、行列の数が多くなればなるほど、あり得る全てのマトリックスの数はどんどん増えていく。それを全て数え上げていくのは、面倒である。そこで、行列の事象が独立という条件のもと、周辺度数を固定した上で、あり得るマトリックスをランダムに出現させる。ただし、全てのマトリックスが出現する確率は全て同じではない。前回の記事でも言及したように、それぞれのマトリックスの出現する確率はフィッシャー大先生によって、公式化されている。よって、ランダムに出やすいマトリックスもあればそうではないマトリックスもある。ただし、もしかしたら1回目の試行で出てきたマトリックスが2回目の試行で出てきたマトリックスと同じになるかもしれないし、そうではないかもしれない。ランダムなので互いの事象は独立である（このあたり、この独立性を緩めたらマルコフ連鎖モンテカルロ法（MCMC）の話になると思われる）。この試行をたくさん繰り返す。Rのfisher.test()のデフォルトは2000回、その試行が行われる。その2000個のマトリクッスのうち、手元にあるデータがどのくらいの確率で出現するのかを計算したのがp値となる。