R:モンテカルロ法で求める円周率のコードを解く。~各務ゼミ生に送る~
おはこんばんにちわ。チャパティです。
この記事は下の記事の続きになります。
ということで、本日も描いていきましょう。
円周率を求めたコード
set.seed(422)
iIter<- 10000000
z <- 0
w <-NULL
pdf( file = "円周率.pdf", height = 5, width = 7, family = "Times")
for( i in 1:iIter)
{
x<-runif(1)
y<-runif(1)
if(x^2+y^2<1.0)
{
z<-z+1
}
if ( i %% 100 == 0)
{
v <- 4.0 * z / i
w <- rbind( w, v )
print( cbind(i,v))
}
}
par( mar = c( 2, 2, 1.5, 0.5))
plot(w,type = "l",xlab = "",ylab = "", bty = "n")
#dev.copy2eps( file = "円周率.eps", height = 5, width = 7, family = "Times")
print(4.0*z/iIter)
dev.off()
これを一行ずつ読み解いていきましょう。
第1ブロック
set.seed(422)
iIter<- 10000000
z <- 0
w <-NULL
pdf( file = "円周率.pdf", height = 5, width = 7, family = "Times")
set.seed(422)
乱数の設定をする。
set.seed(数字)と設定することで毎回同じ乱数を生成することが出来る。
422って書いてるけどここは何でも良し。
iTter<-10000000
z <-0
w <-NULL
それぞれの記号に数字を代入 記号に特に意味はなし
NULLは何もないということ・0ですらない
null を「空欄」的な意味で使う場合があります。 例えば、ゲームで以下のような仕様があったとします。
武器の装備欄は固定で4つある
「1つ目と3つ目に武器を持っていて、2つ目と4つ目には何も持っていない」みたいに、歯抜けがあって、かつ、何番目かに意味がある
https://ufcpp.net/study/csharp/rm_nullusage.html
pdf( file = "円周率.pdf", height = 5, width = 7, family = "Times")
pdfで出力する(詳しくは後ほど)
第2ブロック
for( i in 1:iIter)
{
x<-runif(1)
y<-runif(1)
if(x^2+y^2<1.0)
{
z<-z+1
}
if ( i %% 100 == 0)
{
v <- 4.0 * z / i
w <- rbind( w, v )
print( cbind(i,v))
}
}
for( i in 1:iIter)
{
}
for文
1からiIter(ふつうはnだが、今回は計算式にnが出てきているので便器上iIter)まで繰り返し
要するにiIter回この後の計算をする
x<-runif(1)
y<-runif(1)
runif→一様乱数
xとyに乱数を一つ代入
if(x^2+y^2<1.0)
{
z<-z+1
}
if文
x²+y²<1の時はzにz+1を代入
if ( i %% 100 == 0)
{
}
if文
iを100で割った余りが0の時、{}内のプログラムを実行する
v <- 4.0 * z / i
vに 4.0*z/ i を代入
zはx²+y²≧1の時は0 x²+y²<1の時は1
従って
vはx²+y²≧1の時は0 x²+y²<1の時は4/ i
ここまでで、とりあえず計算終了
円周率の結果は出てくる
w <- rbind( w, v )
rbind→ベクトルを縦に結合する
bindに関して詳しくは下の記事
wのベクトルとvのベクトルを縦に結合したものをwに代入する(plotするときに登場)(がw=NULLなのでどうなるのかよーわからん)(結果的にはx軸がiIterでy軸がvのグラフが出力されてる)
print( cbind(i,v))
i のベクトルとvのベクトルを横に結合したものを出力する
プログラムを走らせたときに
i=100 v=3.2とかが結果として出てくるのはこれを書いてるから
第3ブロック
par( mar = c( 2, 2, 1.5, 0.5))
plot(w,type = "l",xlab = "",ylab = "", bty = "n")
#dev.copy2eps( file = "円周率.eps", height = 5, width = 7, family = "Times")
print(4.0*z/iIter)
dev.off()
par( mar = c( 2, 2, 1.5, 0.5))
par()→図のパラメータを指定
mar = c( 2, 2, 1.5, 0.5)→底辺,左側,上側,右側の順に余白の大きさ
plot(w,type = "l",xlab = "",ylab = "", bty = "n")
plot関数
plot(関数a)→関数aを出力する この場合はw
type = "l"→折れ線で出力する
xlab = "" →x軸の名前を""にする(よーは空白にする)
ylab = "" →y軸の名前を""にする
bty = "n"→グラフの外側の枠を表示しない
print(4.0*z/iIter)
4.0*z/iIterを出力する
print( cbind(i,v))で結果は表示されているものの最後の結果を見るためにこのコードがある
#dev.copy2eps( file = "円周率.eps", height = 5, width = 7, family = "Times")
pdf( file = "円周率.pdf", height = 5, width = 7, family = "Times")
dev.off()
epsでプロットしたものを保存したい場合上のコードを、pdfでプロットしたものを保存したい場合下のコードを使う。
epsの時は最後に1行入れる。pdfの時はif文などで計算を始める前に1行目を入れ、最後にdev.off()で終わらせる。
height→縦の長さ
width→横の長さ
family→文字の書体
まとめ
疲れたー
まだまだ理解してないとこも多いので、これから学んでいかなきゃなあと
次は中心極限定理の証明のプログラムとその解説をすると思います。
それでは、また!
最後まで読んでくれてありがとう!読み終わって内容が面白ければ、「お疲れ様」の意味を込めて「缶コーヒー1杯飲める」程度のサポートをぜひ!