モンテカルロ法による円周率の近似


本屋に情報処理言語のR言語の本がたくさんあって立ち読みしてみたらベクトル処理言語とか書いてあって、あれJ言語みたいじゃんと思いました。
ウィキペディアR言語での「モンテカルロ法による円周率の近似」という例が載っていました。

s <- 100000
x <- runif(s)
y <- runif(s)
sum( x^2+y^2 <= 1)*4/s

記事引用: 『<-』は代入、『runif(a)』は一様乱数をa個作りベクトルで返す関数、『a^2』はaの二乗、『sum(a<=b)』はa<=bであるようなベクトル要素の個数、を意味する。 「サンプル数十万個として、一様乱数でxyのサンプルをつくり、半径1の円弧内に入ったサンプルだけを数える」という計算の本質を反映した記述ができる事が見て取れる。このままRに入力すれば計算が実行できる。以下のように書き換えても意味は同じ。

さて、そのままJ言語に翻訳してみると、次のようになるかと思います。

   s =: 100000
   x =: ? s $ 0
   y =: ? s $ 0
   ( ( +/ ( (x^2) + (y^2)) <: 1) * 4) % s
3.1408

若干というかすごくひいき目だと思いますが、J言語の方がすてきですね!

一応解説ですが、?は乱数発生です。$はシェイプを作りますのでs $ 0でゼロが十万個並んだベクトルができてそれに?を作用させると0~1の乱数が十万個できる。R言語からの翻訳なのでxやyの変数に一度入れているけど、J的には変数使わないと思う。

   ( (+/( ( (? 100000 $ 0)^2)+( (? 100000 $ 0)^2))<:1)*4)%100000
3.14836

しかしこれでは読みにくいですね。あとで考えます。


追記:

   100000 %~ 4 * +/ 1 >: +/ *: ? 2 100000 $ 0
3.1436
   100000 %~ 4 * +/ 1 >: +/ *: ? 2 100000 $ 0
3.13656
   100000 %~ 4 * +/ 1 >: +/ *: ? 2 100000 $ 0
3.14668

解説:かっこがなくなりました! それは見ればわかるか。右から解説すると$(ドルマーク、シェイプ)で2行10万列の0の配列をつくり、それに?(クエスチョンマーク、ランダム)を作用させて20万個の乱数を作ります。*:(アスタリスクコロン)で二乗して、+/(プラススラッシュ)でまず縦にたして10万項目のベクトルになる。>:(だいなりコロン、大なりイコール)を二項動詞として使って1より小さいのを選び+/(プラススラッシュ)で個数を数えます。4倍して、10万で割って終了。%~(パーセントチルダ)は%(パーセント)だけだと左側引数を右側引数で割ることになるので、その順序を~(チルダ)で逆転しています。


こうして考えると、数式を数式として書くのにはJ言語は適していないけれど、試行錯誤しながら効果的に結果を求めるには秀逸な言語ですね。数式を書くのは数式があるので、プログラム言語は数式ではないのだから、J言語のようであるべきだと思うのですが、もう一歩人気ないなぁ。