sicp3.1.2蒙特卡罗方法

#lang sicp
;余数
(remainder 9 7)

;求最大公约数
(define (gcd a b)
  (if (= b 0)
      a
      (gcd b (remainder a b))))

; 线性同余法,a 和 m 是素数
(define (rand-update x)
  (let ((a 48271) (b 19851020) (m 2147483647))
    (modulo (+ (* a x) b) m)))

(define random-init 7)

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

;(rand)
;(rand)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;6/(PI^2) = 两个整数最大公约数为1的概率
;也就是6/概率 = pi^2
;所以pi = (sqrt (/ 6 概率))
;trial试验,试验的次数
;sqrt求平方根
(define (estimate-pi trials)
  (sqrt (/ 6 (monte-carlo trials cesaro-test))))


;返回值为true或者是false,判断两个随机数的最大公约数是否为1(互质)
(define (cesaro-test)
  (= (gcd (rand) (rand)) 1))

;(cesaro-test)
;experiment 试验,最大公约数是否为1
;蒙特卡罗的参数,试验次数,每次试验的值(即:最大公约数是否为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))

;;;;;;;;;;;;;;;;;
(define (estimate-pi-2 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))

;;;;;;;;;;;;;;;;;
(estimate-pi 1000000)
(estimate-pi-2 1000000)




;练习3.5
(define (random-in-range low high)
  (let ((range (- high low)))
    (+ low (random range))))
;    (+ low (* (random range))))) 

(random-in-range 10 20)

(define (estimate-integral P x1 x2 y1 y2 trials)
  (define (experiment)
    (P (random-in-range x1 x2) (random-in-range y1 y2)))
  ;矩形面积
  (let ((rt-area (* (abs (- x1 x2))
                    (abs (- y1 y2)))))
    (* rt-area (monte-carlo trials experiment)))); 矩形面积乘以比率
    
  

THE END
分享
二维码

)">
< <上一篇
下一篇>>