<< 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
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.
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))))
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.
I don't like Shawn solution, I think a predicate should not concern about generating random numbers, that is the job of the experiment. The circle predicate should check whether a point is inside a circle or not, the experiment determines if a random point is inside a predicate, this way you can reuse other predicates.
To martin256. If you look closely, you'll know in Shawn solution, the circle predicate IS the experiment. Just same with the `cesaro-test` in the book.
(define (cesaro-test)
(= (gcd (rand) (rand)) 1))
Abstraction layers can be used here to build more modulized solution making it more suitable to extend to larger system as we have learnt from Chap. 2 Data Abstraction.
2 points to be notified
1. For estimate-integral, instead of hardcode boundaries in test function P, it is better to define P as function taking boundaries as input. Since boundaries define the working region/area to calculate the integral, the concept of working region/area is generic for all kinds of tests in terms of integral estimation. Therefore, it is better to take directly the square/area/region object as parameter instead of x y coordinates. With the objects methods to calculate area and selectors to get boundary points, it is possible to manipulate things at a higher level and thus more intuitive and modulized.
2. For p-estimate, user only needs to input number of trials, as this is the only thing that change the accuracy of the result. Better to hide other paramters by creating abstraction layers with other procedures as we learnt from Chap. 1 Procedural Abstraction.
;; generator and selector for point (define (make-point x y) (cons x y)) (define (xp point) (car point)) (define (yp point) (cdr point)) (define (mid-point p1 p2) (let ((x1 (xp p1)) (x2 (xp p2)) (y1 (yp p1)) (y2 (yp p2))) (define new-x (/ (+ x2 x1) 2)) (define new-y (/ (+ y2 y1) 2)) (make-point new-x new-y) )) ;; generator and selector for circle (define (make-circle pc r) (list pc r)) (define (get-c circ) (car circ)) (define (get-r circ) (cadr circ)) ;; generator and selector for square (define (make-square p1 p2) (list p1 p2)) (define (p1st-square square) (car square)) (define (p2nd-square square) (cadr square)) (define (area-square square) (let ((p1 (p1st-square square)) (p2 (p2nd-square square))) (define x1 (xp p1)) (define x2 (xp p2)) (define y1 (yp p1)) (define y2 (yp p2)) (abs (* (- x2 x1) (- y2 y1))) )) ;; procedure to pick one random point inside square (define rand-in-square (lambda (square) (let ((x1 (xp (p1st-square square))) (x2 (xp (p2nd-square square))) (y1 (yp (p1st-square square))) (y2 (yp (p2nd-square square)))) (make-point (random-in-range x1 x2) (random-in-range y1 y2))))) ;; estimate-integral taking square as input instead of x y coordinates (define (estimate-integral P square trials) (monte-carlo trials (P square))) ;; define unit circle and square used for estimation (define p01 (make-point -1. -1.)) (define p02 (make-point 1. 1.)) (define pm (mid-point p01 p02)) (define unit-square (make-square p01 p02)) (define unit-cir (make-circle pm 1.)) ; helper function to calculate square (define (sqr x) (* x x)) ;; generic in-circle test for point in any square and any circle (define (in-circle? square cir) (lambda () (let ((xc (xp (get-c cir))) (yc (yp (get-c cir))) (r (get-r cir)) (p_rand (rand-in-square square))) (define x_rand (xp p_rand)) (define y_rand (yp p_rand)) (<= (+ (sqr (- x_rand xc)) (sqr (- y_rand yc))) (sqr r))) )) ;; specific test function with unit circle for pi estimate (define in-unit-cir? (lambda (square) (in-circle? square unit-cir))) ;; p-estimate: user only needs to input number of trials (define (p-estimate trials) (* (area-square unit-square) (estimate-integral in-unit-cir? unit-square trials)))
Test:
(p-estimate 1000000) ;;-> 3.141292
---
If you just want to estimate PI, the actual values do not matter. Just take a radius and try to estimate. Actually, which radius you take does also not matter that much (when working with floats)
(define (square x) (* x x)) ; (r * 2)^2 * integral = r^2 * PI ; (2)^2 * integral = PI ; PI = 4 * integral (define (calculatePI integral) (* integral 4.0)) (define (monte-carlo trials experiment) (define (iter trails-remaining trials-passed) (cond ((= trails-remaining 0) (/ trials-passed trials)) ((experiment) (iter (- trails-remaining 1) (+ trials-passed 1))) (else (iter (- trails-remaining 1) trials-passed)))) (iter trials 0)) (define (estimate-integral radius trials) (define squareRadius (square radius)) (define diameter (* 2.0 radius)) (define (predicate) (let ((randomX (random diameter)) (randomY (random diameter))) (<= (+ (square (- randomX radius)) (square (- randomY radius))) squareRadius))) (monte-carlo trials predicate)) (define (estimate-pi radius trials) (calculatePI (estimate-integral radius trials))) (estimate-pi 1.0 100000); 3.13848
---
I tried to make it a little more general.
test runs
(estimate-integral in-circle? 2.0 8.0 4.0 10.0 1000)
=> 3.24
(estimate-integral in-circle? -1.0 1.0 -1.0 1.0 1000)
=> 28.44
(define (in-circle? x1 x2 y1 y2) (let ((center (cons (/ (+ x1 x2) 2.0) (/ (+ y1 y2) 2.0))) (r (/ (min (- x2 x1) (- y2 y1)) 2.0)) (x (random-in-range x1 x2)) (y (random-in-range y1 y2))) (<= (+ (square (- x (car center))) (square (- y (cdr center)))) (square r)))) (define (estimate-integral predicate x1 x2 y1 y2 trials) (let ((area (* (- x2 x1) (- y2 y1)))) (* area (monte-carlo trials (lambda () (predicate x1 x2 y1 y2))))))
---
A specialised solution, not general. Using the random function from the sicp lang in Racket. Just doing 10000 simulations, since this wasn't specified in the problem.
#lang sicp (define (square x) (* x x)) (define (in-unit-circ? x y) (<= (+ (square x) (square y)) 1)) (define (rand-between x1 x2) (+ x1 (random (- x2 x1)))) (define (estimate-integral P x1 x2 y1 y2) ; x1 < x2 && y1 < y2 (define area (abs (* (- x2 x1) (- y2 y1)))) (* area (monte-carlo 10000 (lambda () (P (rand-between x1 x2) (rand-between y1 y2)))))) ;; (estimate-integral in-unit-circ? -1.0 1.0 -1.0 1.0) ;; => 3.148
Using the monte-carlo procedure as defined in the book, and the exact->inexact primitive procedure to turn a fraction into a float: