乐正

Actions speak louder than words.

Sicp-ex3-5

问题

蒙特卡洛积分是一种通过蒙特卡洛模拟估计定积分值的方法。考虑由谓词$P(x, y)$描述的一个区域的面积计算问题,该谓词对于此区域内部的点$(x, y)$为真,对于不在区域内的点为假。举例来说,包含在以$(5, 7)$为圆心半径为3的圆圈所围成的区域,可以用检查公式$(x-5)^2 + (y-7)^2 <= 3^2$是否成立的谓词描述。要估计这样一个谓词所描述的区域的面积,我们应首先选取一个包含该区域的矩形。例如,以$(2, 4)$和$(8, 10)$为对角线的矩形包含着上面的圆。需要确定的积分也就是这一矩形中位于所关注区域内的那个部分。我们可以这样估计积分值:随机选取位于矩形中的点$(x, y)$,对每个点检查$P(x, y)$,确定该点是否位于所考虑的区域内。如果试了足够多的点,那么落在区域内的点的比率能够给出矩形中有关区域的比率。这样,用这一比率去乘整个矩形的面积,就能得到相应积分的一个估计值。

将蒙特卡洛积分实现为一个过程estimate-integral,它以一个谓词$P$,矩形的上下边界$x1$、$x2$、$y1$和$y2$,以及为产生估计值而要求试验次数作为参数。你的过程应该使用上面用于估计$\pi$值的同一个monte-carlo过程。请用你的estimate-integral,通过对单位圆面积的度量产生出$\pi$的一个估计值。

你可能发现,有一个从给定区域中选取随机数的过程非常有用。下面的random-in-range过程利用1.2.6节中的random实现这一工作,它返回一个小于其输入的非负数。

1
2
3
(define (random-in-range low high)
  (let ((range (- high low)))
    (+ low (random range))))

解答

假设单位圆是圆心在$(1, 1)$的单位圆:

练习3.5 (ex3-5.scm) download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
(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 (random-in-range low high)
  (let ((range (- high low)))
    (+ low (random range))))

(define (estimate-integral p x1 x2 y1 y2 trials)
  (define (test)
    (p (random-in-range x1 y2)
       (random-in-range x1 y2)))
  (* (- x2 x1) (- y2 y1) (monte-carlo trials test)))

(define (in-range? x y)
  (<= (+ (square (- x 1))
         (square (- y 1)))
      1))

测试

1
2
3
4
5
6
7
8
(estimate-integral in-range? 0 2.0 0 2.0 1000)
;Value: 3.112

(estimate-integral in-range? 0 2.0 0 2.0 10000)
;Value: 3.1404

(estimate-integral in-range? 0 2.0 0 2.0 100000)
;Value: 3.14968

draft

« sicp-ex3-4 sicp-ex3-6 »

Comments