モンテカルロ法で定積分して面積求めるとかやります.
汎用なモンテカルロ法はすでに用意されてて,こんなのです.
trialsが試行回数で,experimentは実行する関数(predicate)だ!で,このmonte-carloはexperimentがどれだけの割合でtrueを返したかの割合を返してくれる.
(define (monte-carlo trials experiment) (define (iter trials-remaining trials-passed) (cond ((= trials-remaining 0) (/ trials-passed trials)) ((experiment) (iter (- trials-remaining 1) (+ trials-passed 1))) (else (iter (- trials-remaining 1) trials-passed)))) (iter trials 0))
んじゃ定積分はどうするのかというと,まず積分領域を内包する長方形を作る.
そしてこの長方形上にあるランダムな点(x,y)を決めます.
で,その点が積分領域上にあればtrue,そうでなければfalseにするpredicateをでっちあげて,積分領域上に点があるかを判断する.
これをムチャクチャ繰り返すと,どれだけの割合で積分領域上にあったのか,その割合が出る.その割合に長方形の面積をかければ,それが定積分の近似値です!
[0,n]の乱数はgaucheだと,こんな感じで作れます.
(use srfi-27) (define (random n) (* (random-real) n))
これを利用して,[low, high]内で一様に分布した乱数はrandom-in-rangeとしてこんな風に定義する.
(define (random-in-range low high) (let ((range (- high low))) (+ low (random range)))) (random-in-range 1 5) ; 3.35598268950961 (random-in-range 1 5) ; 4.653520179856983 (random-in-range 1 5) ; 3.415157581535672
近似値を求める関数はこれ.
(x1, y1),(x2, y2)を結ぶ線分が対角線になるような長方形が積分領域を内包する長方形になる.
Pは2つの引数(x,y)をとって,それが積分領域にあるかどうかを判断するpredicateで,trialsが試行回数です.
(define (estimate-integral P x1 x2 y1 y2 trials) (let ((low-x (min x1 x2)) (low-y (min y1 y2)) (high-x (max x1 x2)) (high-y (max y1 y2))) ; 長方形内のランダムな点をとって,それが積分領域内かどうかを返す (define (exp-point-in-region) (let ((rand-x (random-in-range low-x high-x)) (rand-y (random-in-range low-y high-y))) (P rand-x rand-y))) ; 主部 (* (monte-carlo trials exp-point-in-region) (- high-x low-x) (- high-y low-y))))
課題は,を定積分で求めることなので,長方形は(-1, -1),(1,1)を対角線とする長方形,積分領域は(0, 0)を中心とする単位円にします.そういうわけで,predicateはこうなる.
(define (square x) (* x x)) (define (pred-in-unit-circle x y) (< (+ (square x) (square y)) 1))
あとは実行するだけですね.
(estimate-integral pred-in-unit-circle -1 1 -1 1 1000000) ; 3.141024