モンテカルロ法で円周率を求める
モンテカルロ法で円周率を求める
投稿日2017.8.28
対象となるある現象がいくつかの確率変数の関数で表現できる場合（数値モデル）、各変数で仮定される確率分布に沿った標本値を大量に生成（乱数）して、その計算結果の分布を推定（結果の推定）する方法です。
数値モデル → 乱数挿入 → 結果の推定
モンテカルロ法を使用して円周率 (π) を求めてみます
上図の四角の中にランダムに点を打ってみます
四角の面積は 2×2=4、円の面積は π × 1^2 = π
打った点の個数をn、円の中に入った個数をmとします
非常に多くの点を打つことで四角の中は塗りつぶされます
そこで次のような比率を考えます
点を多く打つことで、全体の点の数と円内の点の数との比が
四角の面積と円の面積の比に等しくなるはずです．
n : m = 4 : π
π＝4 × m / n （数値モデル）
Rを使用してnの乱数を多く作成すれば円周率が求められることになります
ではやってみます
一様乱数から(x,y)を作成します （乱数）
-1<x<1, -1<y<1
試しに5組分だけ作成してみます
dat<-matrix(runif(5*2,-1,1),2)
イメージ
一行目をｘ、二行目をｙとして x^2 + y^2 < 1 となる点の個数が
円内の個数mということになります
m<-colSums(dat^2) < 1
全ての点をプロットしてｍは赤で示します
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
m はcolSums(dat^2) < 1の場合にはTRUE＝1、それいがいはFALSE＝0 ここがポイント
col=(2-m)とした場合には円内に入ったら1（黒）、円外では2（赤）となります
これまでの式をまとめて100個打ってみます
dat<-matrix(runif(100*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
分かりにくいので10000個打ってみます
dat<-matrix(runif(10000*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
おっ！円が見えてきました！！！
ちょっと感動
dat<-matrix(runif(n*2,-1,1),2)
m<-colSums(dat^2) < 1
plot( dat [1,], dat [2,], col = (2-m) , pch = 20 )
結果の推定
π＝4 × m / n を推定します
上記の結果より
π＝4 × m / 10000 = 3.1388 (乱数なので結果は異なります)
まとめて関数にしてみます
pai<-function(n){
dat<-matrix ( runif ( n*2 , -1 , 1 ) , 2 )
m<-colSums( dat^2 ) < 1
pl<-plot( dat [1,] , dat [2,] , col = (2-m) , pch = 20 , xlab = n , ylab = "" )
return( ( 4*sum( m ) ) / n )
}
nを100,500,1000,5000と変更して図を描いてみます
par(mfrow = c(2,2))
pai(100) ;pai(500) ;pai(1000) ;pai(5000)
dev.off()
n=100 π＝ 3.4
n=500 π＝ 3.128
n=1000 π＝ 3.144
n=5000 π＝ 3.1352
二元配置分散分析（各水準の信頼区間推定）
二元配置分散分析（各水準の信頼区間推定）
投稿日2017.8.24
各要因における各水準の母平均の区間推定
繰り返し数 n
n = 繰り返し数
要因1水準Aの有意水準αの信頼区間は、xai ± t α/2 (自由度) * 標準誤差 となります． 母分散が未知なので、標準誤差を残差の平均平方（不偏分散） / 水準のサイズ としてます．したがってｔの自由度は、残差の平均平方 の自由度になります．
95％信頼区間
xai ± t 0.025 ( f*g*(n-1) ) * 標準誤差
例）次の水準Aの95％信頼区間を求めてみます
xai
= mean(c(s1,s2,s3))
= 13.1333
自由度2*3*(5-1)=24, 有意水準5％のｔ値は t(24)(0.025)
qt ( 0.975 ,24 ,lower.tail=TRUE )
= 2.063899
残差の平均平方（不偏分散）は3.08
残差の標準誤差は√(3.08/5*3) = √(3.08/15)
したがってxa1.の95%信頼区間は
mean ( c (s1,s2,s3) ) - 2.063899 * sqrt ( 3.08/15 )
= 12.1981
mean ( c (s1,s2,s3) ) + 2.063899 * sqrt ( 3.08/15 )
= 14.06856
次のような式で表現できます
交互作用がない場合
Ai水準での母平均の推定量は xai.
Ai水準での95％信頼区間は xai. ± t 0.025 ( de(e) ) × √ ( V(e) / g×n )
Bi水準での母平均の推定量は xa.j
Ai水準での95％信頼区間は xa.j ± t 0.025 ( de(e) ) × √ ( V(e) / f×n )
交互作用がある場合
(AiBj) での母平均の推定量は xaij
(AiBj) での95％信頼区間は xaij ± t 0.025 ( de(e) ) × √ ( V(e) / n )
トービット回帰直線
トービット回帰直線
投稿日2017.8.4
打ち切りデータの場合、つまり天井効果や床効果が生じているデータの場合には、トービット回帰直線で分析します．
参考書はもちろん
豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012
Rのサンプルanscombeを使います．( n=11の4セット )
x1<-anscombe$x1
x4<-anscombe$x4
y1<-anscombe$y1
y4<-anscombe$y4
par(mfrow = c(1,2))
plot(x1,y1,xlab = "図A",ylab = "")
plot(x4,y4,xlab = "図B",ylab = "")
dev.off()
図A 通常の散布図
lm1<-lm(y1~ x1,data= anscombe)
plot(x1,y1)
abline(lm1$coefficients[1],lm1$coefficients[2])
視覚的にも何となくこの直線でフィットしているように見えます
ｘ軸を予測値にして、残差を標準化して更に見やすくしてみます
ei<- y1-( lm1$coefficients[1]+ lm1$coefficients[2]* x1)
eis<-ei/sqrt( (sum(ei^2) )/(11 - 2) ) #標準化残差
plot(lm1$coefficients[1]+ lm1$coefficients[2]* x1, eis,xlab="図A（ｘ軸：予測値）",ylab=" ", xaxt="n")
abline(h=0, lty=2, col=2)
図B 打ち切りデータがある場合
トービット回帰直線（ 関数vglmを使用 ）
通常の回帰直線（実践）
y4x4 <- lm(y4 ~ x4); summary(y4x4)
plot(x4,y4)
abline(summary(y4x4)$coefficients1,summary(y4x4)$coefficients2)
関数vglmを使用してトービット回帰直線を破線で加えます
y4x4t <- vglm(y4 ~ x4, tobit(Lower=8, Upper=19), trace=TRUE); summary(y4x4t)
abline(summary(y4x4t)@coefficients1,summary(y4x4t)@coefficients3,lty=2)
床効果で下に引っ張られていた回帰直線をやや上方修正できました．
