SICP 3.1.2 The Benefits of Introducing Assignment
章の最初は軽くて快調、快調。と思いきや乱数の「質」は意外と重要である事を再確認するはめに。
ノート
代入を導入する事によるメリット。
As we shall see, introducing assignment into our programming language leads us into a thicket of difficult conceptual issues. Nevertheless, viewing systems as collections of objects with local state is a powerful technique for maintaining a modular design.
代入（と言うか再代入）が概念的に難しい問題となる。それでもシステムを状態変数を持ったオブジェクトの集まりと捉える事が部品化された設計を維持する強力な手段である。
前の値を引数に取り次の乱数を返すrand-updateがあるとして、この「前回値」を保存して毎回rand-updateに渡す機能をrandが抱え込む事でrand-update固有の作業をアプリケーション側にさせず、乱数発生とモンテカルロ法の部分を分離する事が出来る。
脚注６にはrand-updateの実装法のひとつとして線形合同法が挙げられているが、実際にやってみるとなかなか良い結果は得られない。英語版WikipediaのLinear congruential generatorの項に各種の定数の例が出ているが、どれを試しても得られるπの値は2.7近辺。たまたま見つけたa=2100005341、b=1、m=4294967087だけがそれなりの結果を出した。
(define random-init 103) (define (rand-update x) (let ((a 2100005341) (b 1) (m 4294967087)) (remainder (+ (* a x) b) m)))
rand-updateを直接呼ぶ方はestimate-pi2と名前を変えて動作確認。
> (estimate-pi 10000) 3.1452865847037783 > (estimate-pi2 10000) 3.1452865847037783 >
初期値が同じなので同じ答えになる。
そもそもWikipedia（英語版）には線形合同法はモンテカルロ法には適さないと書いてある。
Raketのrandomを使えばかなり良い結果が得られる。
(define random-init 1) (define (rand-update x) (random 4294967087))
動かすと
> (estimate-pi 10000) 3.1367645067819305 > (estimate-pi2 10000) 3.135222467603063 >
こちらは毎回異なる数列となる為、実行する度に答えが変わる。
Exercise 3.5
まずrandom-in-rangeで呼び出すrandomがこのままでは整数しか返さない。
> (random-in-range -1 1) 0 > (random-in-range -1 1) -1 > (random-in-range -1 1) -1 > (random-in-range -1 1) 0 > (random-in-range -1 1) 0 > (random-in-range -1 1) -1 >
なのでRacketのrandomを使って以下の様に変更する。
(define (random-in-range low high) (let ((range (- high low))) (+ low (* (random) range))))
これで以下の様な動作になる。
> (random-in-range -1 1) 0.9767107966253177 > (random-in-range -1 1) -0.9212087592136632 > (random-in-range -1 1) 0.4519761470172179 > (random-in-range -1 1) 0.5802750309689919 > (random-in-range -1 1) -0.12653579011546545 > (random-in-range -1 1) 0.25275783673246166 > (random-in-range -1 1) 0.735915879502544 >
Pとは別にx,yの範囲が渡されると言う事はP自身はrandom-in-rangeを呼ばない様にしないとコードが冗長になってしまう。なのでPには任意の点(x, y)を引数として渡す事にする。monte-carloはtrialに対してtrueとなった「割合」が返って来るので、(x1, y1)〜(x2, y2)の面積との積を取る。ここで答えが分数にならない様に1.0を掛けておく。
(define (estimate-integral P x1 x2 y1 y2 trials) (define (test) (P (random-in-range x1 x2) (random-in-range y1 y2))) (* (* (- x2 x1) (- y2 y1)) (monte-carlo trials test) 1.0))
πは半径１の面積。なので単位円の内側か否かをテストする。
(define (in-unit-circle x y) (< (sqrt (+ (* x x) (* y y))) 1.0)) (define (estimate-pi trials) (estimate-integral in-unit-circle -1 1 -1 1 trials))
動作確認
> (estimate-pi 100) 3.2 > (estimate-pi 1000) 3.112 > (estimate-pi 10000) 3.1404 > (estimate-pi 100000) 3.1444 > (estimate-pi 1000000) 3.141944 >
Exercise 3.6
- 各呼び出しで共通のxを使う為にはクロージャを返す必要がある。
- そのクロージャはメッセージとしてアトムを受け取る。
- 'generateを受け取ったらxを乱数で入れ替えてxを返す。
- 'resetを受け取ったら引数をひとつ取り、その引数でxを入れ替える関数を返す。
(define rand (let ((x random-init)) (lambda (msg) (cond ((eq? msg 'generate) (begin (set! x (rand-update x)) x)) ((eq? msg 'reset) (lambda (val) (set! x val))) (else (error "unknown message"))))))
実行する。
> (rand 'generate) 1552195774 > ((rand 'reset) 100) > (rand 'generate) 3842113925 > ((rand 'reset) 103) > (rand 'generate) 1552195774 > (rand 'dummy) . . unknown message >
シードを元に戻せば最初の乱数列に戻る。
