sicp-ex-3.5


<< Previous exercise (3.4) | Index | Next exercise (3.6) >>


Given:

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

The solution is:

 (define (P x y) 
   (< (+ (expt (- x 5) 2) 
         (expt (- y 7) 2)) 
      (expt 3 2))) 
 (define (estimate-integral P x1 x2 y1 y2 trials) 
   (define (experiment) 
     (P (random-in-range x1 x2) 
        (random-in-range y1 y2))) 
   (monte-carlo trials experiment)) 

Test:

 (estimate-integral P 2.0 8.0 4.0 10.0 100) 

Then we can estimate pi with the fact that a circle area is (pi * r²).

Hence pi ≅ (Monte Carlo results * rectangle area) / r²

 (define pi-approx 
   (/ (* (estimate-integral P 2.0 8.0 4.0 10.0 10000) 36) 
      9.0)) 
 pi-approx 

Which gave 3.1336 during my test.

This function has to be tested under MIT Scheme, neither gambit-scheme or SISC implements (random) - actually (random) is not part of R5RS nor SRFI.

NB: using 2.0 instead of 2 in estimate-integral is primordial. If you pass two integers to (random-in-range low high), it will return another integer strictly inferior to your 'high' value — and this completely screws the Monte-Carlo method (it then estimates pi to ~3.00).


Using Racket's built-in (random), which returns a float between 0 and 1, the following set of functions can be used to find the numerical integral of a function with a predicate, using the provided (monte-carlo) function:

 (define (random-in-range low high) 
   (let ([range (- high low)]) 
     (+ low (* (random) range)))) 
  
 (define (get-area lower-x lower-y upper-x upper-y) 
   (* (- upper-x lower-x) 
      (- upper-y lower-y))) 
  
 (define (estimate-integral pred lower-x lower-y upper-x upper-y trials) 
   (let ([area (get-area lower-x lower-y upper-x upper-y)] 
         [experiment (lambda () (pred (random-in-range lower-x upper-x) 
                                      (random-in-range lower-y upper-y)))]) 
     (* area (monte-carlo trials experiment)))) 

Testing with a large number of trials:

 (estimate-integral (lambda (x y) (<= (+ (* x x) (* y y)) 1.0)) -1.0 -1.0 1.0 1.0 100000000) 
 3.14130636 

Shawn

I think the procedure estimate-integral is not very nicely formulated in the text above, in my opinion, a better solution for this exercise would be:

 (define (estimate-integral P x1 x2 y1 y2 trials) 
   (* (* (- x2 x1) 
         (- y2 y1)) 
      (monte-carlo trials P))) 
  
 (define (in-circle) 
   (>= 1 (+ (square (random-in-range -1.0 1.0)) 
            (square (random-in-range -1.0 1.0))))) 
  
 (define (random-in-range low high) 
   (let ((range (- high low))) 
     (+ low (random range)))) 
  
 (define (estimate-pi) 
   (estimate-integral in-circle -1.0 1.0 -1.0 1.0 1000)) 
  
 ;; monte carlo procedure 
 (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)) 

We should leave the job of generating random numbers to the predicate, instead of intertwining it with estimate-integral. If we do it like this, we cannot adapt estimate-integral to other predicates.


athird

If you're using Racket the random function doesn't work as described in the text. You might want to use this instead:

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

gambiteer

Gambit has (random-real) for a random inexact real between 0 and 1 and (random-integer n) for a random integer between 0 and n.