3.1.2
(define (estimate-pi trials)
(sqrt (/ 6 (monte-carlo trials cesaro-test))))
(define (cesaro-test)
(= (gcd (rand) (rand)) 1))
(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)
)
; randのかわりにrand-updateを使うようにする
(define (estimate-pi trials)
(sqrt (/ 6 (random-gcd-test trials random-init))))
(define (random-gcd-test trials initial-x)
(define (iter trials-remaining trials-passed x)
(let ((x1 (rand-update x)))
(let ((x2 (rand-update x1)))
(cond
((= trials-remaining 0) (/ trials-passed trials))
((= (gcd x1 x2) 1)
(iter
(- trials-remaining 1)
(+ trials-passed 1)
x2
)
)
(else
(iter (- trials-remaining 1) trials-passed x2))))))
(iter trials 0 initial-x)
)
gaucheにはrand-initないからこのままでは動かない。
Ex 3.5
モンテカルロ積分を、⼿続き estimate-integral として実装せよ。
引数として、述語 P、⻑⽅形の上下界 x1, x2, y1, y2、推定値を⽣成するための試⾏回数を取る。
この⼿続きは、上で π を推定するために使ったものと同じ monte-carlo ⼿続きを使わなければならない。
この estimate-integral を使って、単位円の⾯積を測ることによって π の推定値を求めよ。
(use srfi-27)
; このモジュールはメルセンヌツイスタアルゴリズム (Mersenne Twister乱数発生器 参照) を基礎に用いた
; SRFI-27疑似乱数発生器インタフェースを提供します。
(print "==ex3.5==")
(define (random-in-range low high)
(let ((range (- high low)))
(+ low (* range (random-real))))) ; gaucheのrandom-integerを使う
(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.0)
)
(define (estimate-integral P x1 x2 y1 y2 trials)
(*
(* (- x2 x1) (- y2 y1))
(monte-carlo trials (lambda () (P (random-in-range x1 x2) (random-in-range y1 y2)))))
)
;テスト手続き
(define (ex-3.5)
;; 中心(5, 7) 半径3 の円の場合
(define (p-test x y)
(<
(+ (expt (- x 5) 2) (expt (- y 7) 2))
(expt 3 2)))
;; (0, 0)の半径1の円 (piの推定)
(define (pi-test x y)
(<
(+ (expt x 2) (expt y 2))
(expt 1 2)))
(print "==== answer " (* 3.1415 (* 3 3)))
(print "100回: " (estimate-integral p-test 2 8 4 10 100))
(print "10000回: " (estimate-integral p-test 2 8 4 10 10000))
(print "100000回: " (estimate-integral p-test 2 8 4 10 100000))
(print "==== answer " 3.1415)
(print "100回: " (estimate-integral pi-test -1 1 -1 1 100))
(print "10000回: " (estimate-integral pi-test -1 1 -1 1 10000))
(print "100000回: " (estimate-integral pi-test -1 1 -1 1 100000))
(print "1000000回: " (estimate-integral pi-test -1 1 -1 1 1000000))
;==ex3.5==
;==== answer 28.273500000000002
;100回: 27.36
;10000回: 28.2852
;100000回: 28.26144
;==== answer 3.1415
;100回: 2.76
;10000回: 3.1144
;100000回: 3.14936
;1000000回: 3.141632
)
(ex-3.5)
Ex 3.6
乱数⽣成器をリセットして、与えられた値から始まる数列を作れるようにできると便利だ。
記号 generate またはreset のどちらかを引数として取り、次のようなふるまいをする新しい rand ⼿続きを設計せよ。
(rand 'generate) は新しい乱数を⽣成する。
((rand 'reset) ⟨new-value ⟩) は内部の状態変数を指定された ⟨new-value⟩ にリセットする。
つまり、状態をリセットすることによって、再現可能な数列が⽣成できる。
これは乱数を使うプログラムをテストしたりデバッグしたりするのにとても役に⽴つ。
=> gaucheのrandom seedの設定（random-source-pseudo-randomize!）は引数を２つとるのでふたつ で実装。
new-valueから始まる・・・にはしてない(´・ω・｀)
;gaucheのランダム
;random-source-pseudo-randomize! s i j
;https://practical-scheme.net/gauche/man/gauche-refj/randamubitutonososu.html
(use srfi-27)
(define (rand arg)
(cond
((eq? arg 'generate) (random-real))
((eq? arg 'reset)
(lambda (i j) (random-source-pseudo-randomize! default-random-source i j)))
(else "Failed"))
)
(define (ex-3.6)
((rand 'reset) 1 2)
(print (rand 'generate))
((rand 'reset) 3 4)
(print (rand 'generate))
((rand 'reset) 1 2)
(print (rand 'generate))
;0.6871964570424902
;0.726033823048644
;0.6871964570424902
)
(ex-3.6)