sicp3.1.2蒙特卡罗方法_monte carlo trials-CSDN博客

这篇博客探讨了如何利用Scheme编程语言实现求余数、最大公约数（gcd）的递归算法，并应用线性同余法生成随机数。通过蒙特卡罗方法，作者估算出了两个随机数的最大公约数为1的概率，进而推算π的值。实验中，通过大量试验来逼近6/π²的比值，从而得到π的近似值。此外，还介绍了如何在程序中嵌入随机数生成以进行更精确的计算。

摘要由CSDN通过智能技术生成

#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)))); 矩形面积乘以比率