sicp 3章 問題3.6 まで。2013-05-13 担当分。

q3.5.scm

(use slib)
(load "/usr/local/slib/random")
(load "./s3.1.2_estimate-pi-using-rand.scm")

(use srfi-27) ; random-real

(print "---------------------------- q3.5 -----------------------------")

; モンテカルロ積分
; たとえば四角形をちょうど半分にし、その片方の面積を求めたい場合、
; ランダムに点を取り、その点が求めたい方に入っているかどうかの確率は1/2
; つまり、(四角形の面積) * 1/2 で求めたい面積が求まる

(define (random-in-range low high)
  (let ((range (- high low)))
    (+ low (* (random-real) range)))) ; (+ low (random range))))

(define (estimate-integral P x1 x2 y1 y2 trials)
  (define (dist a1 a2)
    (abs (- a1 a2)))
  (define (area x1 x2 y1 y2)
    (* (dist x1 x2) (dist y1 y2)))
  (* (area x1 x2 y1 y2)
     (monte-carlo trials
                  (lambda ()
                    (define x (rand-in-rang x1 x2))
                    (define y (rand-in-rang y1 y2))
                    (P x y)))))


(use gauche.test)

(test-start "estimate-integral")

(define (rand-in-rang a1 a2)
  (if (>= a2 a1)
    (random-in-range a1 a2)
    (random-in-range a2 a1)))

(test-section "case: square using corners at (0, 0) and (100, 100)")
(define x1 0)
(define x2 200)
(define y1 0)
(define y2 200)
(define (in-square? x y)
  (define X1 0)
  (define X2 100)
  (define Y1 0)
  (define Y2 100)
  (and (>= x X1) (<= x X2) (>= y Y1) (<= y Y2)))
(test "is erea 10000? (margin:200)" 10000 (lambda () (estimate-integral in-square? x1 x2 y1 y2 1000))
      (lambda (expected actual) (<= (abs (- expected actual)) 200)))

(test-section "case: circle using corners of radius 100 centered (100, 100)")
(define x1 0)
(define x2 200)
(define y1 0)
(define y2 200)
(define (in-circle? x y)
  (define R 100)
  (define X 100)
  (define Y 100)
  (>= (* R R) (+ (* (- x X) (- x X)) (* (- y Y) (- y Y)))))
(test "is erea about 31400? (margin:200)" 31400 (lambda () (estimate-integral in-circle? x1 x2 y1 y2 1000))
      (lambda (expected actual) (<= (abs (- expected actual)) 200)))

(test-end)

(test-start "estimate-pi using estimate-integral")

(define x1 0)
(define x2 2)
(define y1 0)
(define y2 2)
(define (in-unit-circle? x y)
  (define R 1)
  (define X 1)
  (define Y 1)
  (>= (* R R) (+ (* (- x X) (- x X)) (* (- y Y) (- y Y)))))
(test "is erea about 3.14?" 3.14 (lambda () (estimate-integral in-unit-circle? x1 x2 y1 y2 1000))
      (lambda (expected actual) (<= (abs (- expected actual)) 0.2)))

(test-end)

;(print (* (estimate-integral in-unit-circle? x1 x2 y1 y2 20000000) 1.0))

q3.6.scm

; s3.1.2_rand.scmと同じでrand自体はlambdaを返し、その引数がsymbolになっていれば良さそう

(define a 3)
(define b 5)
(define m 13)
(define random-init 8)
(define (rand-update x)
  (modulo (+ (* a x) b) m))

(define rand
  (let ((x random-init))
    (define (rand-generate)
      (set! x (rand-update x))
      x)
    (define (rand-reset)
      (lambda (new-x)
        (set! x new-x)))
    ; lambdaでもいいよ
    ;(define dispatch (lambda (command)
    ;                   (cond ((eq? command 'generate) (rand-generate))
    ;                         ((eq? command 'reset) (rand-reset))
    ;                         (else error "no command"))))
    (define (dispatch command)
      (cond ((eq? command 'generate) (rand-generate))
            ((eq? command 'reset) (rand-reset))
            (else error "no command")))
    dispatch))

(use gauche.test)
(test-start "new rand")

(test-section "case: generate")
(test* "rand?" 3 (rand 'generate))
(test* "rand?" 1 (rand 'generate))
(test* "rand?" 8 (rand 'generate))
(test* "rand?" 3 (rand 'generate))
(test* "rand?" 1 (rand 'generate))
(test* "rand?" 8 (rand 'generate))

(test-section "case: reset")
; 3 -> 1 -> 8 -> 3 -> 1 -> 8
(rand 'generate) ;=> 3
((rand 'reset) 1)
(test* "rand?" 8 (rand 'generate))

(test-end)

s3.1.2_estimate-pi-using-rand-update.scm

; （代入がないとすると）s3.1.2_estimate-pi-using-rand.scm と違い、
; rand-updateを直接使うしかない
; つまり、iter で回す時に乱数を覚えておかなければいけない

(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))

s3.1.2_estimate-pi-using-rand.scm

(load "./s3.1.2_rand.scm")

; ランダムな2つの数の gcd が 1 になる確率は 6/(pi^2)
; なので、pi は逆数とって平方根取ればよい
(define (estimate-pi trials)
  (sqrt (/ 6 (monte-carlo trials cesaro-test))))

(define (cesaro-test)
   (= (gcd (rand) (rand)) 1))

; trials 分繰り返し、最後に (experiment が #t の回数) / traials の確率を出す
(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))


(print "--------------------------------------------------------------------------------------")
(print (estimate-pi 100))

s3.1.2_rand.scm

; refs. http://ja.wikipedia.org/wiki/%E7%B7%9A%E5%BD%A2%E5%90%88%E5%90%8C%E6%B3%95

(define a 3)
(define b 5)
(define m 13)
(define random-init 8)
(define (rand-update x)
  (modulo (+ (* a x) b) m))

(define rand
  (let ((x random-init))
    (lambda ()
      (set! x (rand-update x))
      x)))

(use gauche.test)
(test-start "rand")
(test* "rand?" 3 (rand))
(test* "rand?" 1 (rand))
(test* "rand?" 8 (rand))
(test* "rand?" 3 (rand))
(test* "rand?" 1 (rand))
(test* "rand?" 8 (rand))
(test-end)