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は何もないということ・０ですらない

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文
１から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²＜１の時はzにz+1を代入

if ( i %% 100 == 0)
 {
 
 }

if文
iを100で割った余りが０の時、｛｝内のプログラムを実行する

v <- 4.0 * z / i

ｖに 4.0＊ｚ/ i を代入
ｚはx²＋y²≧1の時は０　x²＋y²＜１の時は1
従って
ｖはx²＋y²≧1の時は０　x²＋y²＜１の時は４/ i 

ここまでで、とりあえず計算終了

円周率の結果は出てくる

w <- rbind( w, v )

rbind→ベクトルを縦に結合する

bindに関して詳しくは下の記事

Rでのデータの編集と演算

ｗのベクトルとｖのベクトルを縦に結合したものをｗに代入する（plotするときに登場）(がw=NULLなのでどうなるのかよーわからん)（結果的にはｘ軸がiIterでy軸がvのグラフが出力されてる）

print( cbind(i,v))

 i のベクトルとｖのベクトルを横に結合したものを出力する

プログラムを走らせたときに

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→文字の書体

まとめ

疲れたー

まだまだ理解してないとこも多いので、これから学んでいかなきゃなあと

次は中心極限定理の証明のプログラムとその解説をすると思います。

それでは、また！