2006/08/27

いまだゴールが見えないSICP。ようやく第3章が終わった。思い返せば第2章を終えたのが今年の3月22日で、それから5ヵ月もかかってる。おそすぎ。まあゆっくりやります。

最後の問題は、また乱数とモンテカルロだった。ところで第3章では疑似乱数を線形合同法で求める方法が紹介されているので(226ページの注6)、Gauche に用意されているメルセンヌ・ツイスターを使っちゃずるいと思う。いや、べつに誰がずるいというわけでなく、それでは自分的に勉強にならないと思う。ってな大口を叩いてたところで、Knuth本を持ってないからWikipediaに載っている係数をひろってくるだけなんだけど。
(define (monte-carlo experiment-stream passed failed)
(define (next passed failed)
(stream-cons
(/ passed (+ passed failed))
(monte-carlo
(stream-cdr experiment-stream) passed failed)))
(if (stream-car experiment-stream)
(next (+ passed 1) failed)
(next passed (+ failed 1))))

;; ex. 3.81
(define rand-init 1)
(define (rand-update x)
(modulo (* x 1664525) 1013904223))
(define (random-numbers c)
(define (rands init)
(stream-cons init
(stream-map rand-update
(stream-delay (rands init)))))
(cond ((equal? c 'generate)
(rands rand-init))
((equal? c 'reset)
(lambda (init)
(rands init)))
(else
(error "Unknown argument" c))))

;; ex. 3.82
(define (random-numbers-in-range x1 x2 rand-init)
(define (normalize x)
(+ (* x (/ (- x2 x1) 1013904223)) x1))
(stream-map normalize ((random-numbers 'reset) rand-init)))
(define (estimate-integral P x1 x2 y1 y2)
(stream-scale
(monte-carlo
(stream-map P
(random-numbers-in-range x1 x2 1)
(random-numbers-in-range y1 y2 506952111))
0 0)
(* (- x2 x1) (- y2 y1))))

(define (P x y)
(<= (+ (square (- x 5))
(square (- y 7)))
9))
(define ex-3-82
(estimate-integral P 2 8 4 10))
(define (ex-3-82-pi howfar)
(/ (stream-ref ex-3-82 howfar) 9))

gosh > (ex-3-82-pi 1000)
3.1128871128871127
最初、x軸とy軸について同じ random-numbers-in-range を使ってたら、どうにも想定している結果よりかなり小さい値しか求まらない。そりゃあ y=x 上の点しか拾ってないんだから、よく考えれば当り前の話だった。結果的に初期値として選んだ 1 と 506952111 は適当で、あまり背景はない(1013904223の約半分)。これくらいオーダーを違くすれば、そこそこバラけるような。

No comments :