CODEWITHSUNDEEP
schememontecarlopi
John Friedrich
John Friedrich

Reputation: 343

Scheme Monte-Carlo-Sampling

I am trying to determine the number of marbles that fall within a given circle (radius 1) given that they have random x and y coordinates.

My overall goal is to find an approximate value for pi by using monte carlo sampling by multiplying by 4 the (number of marbles within the circle)/(total number of marbles).

I intended for my function to count the number of marbles within the circle, but I am having trouble following why it does not work. Any help on following the function here would be appreciated.

Please comment if my above request for help is unclear.

(define(monte-carlo-sampling n)
 (let ((x (- (* 2 (random)) 1))
       (y (- (* 2 (random)) 1)))
 (cond((= 0 n) 
    * 4 (/ monte-carlo-sampling(+ n 1) n) 
     ((> 1 n) 
     (cond((< 1 (sqrt(+ (square x) (square y))) (+ 1 (monte-carlo-sampling(- n 1)))))
         ((> 1 (sqrt(+ (square x) (square y))) (monte-carlo-sampling(- n 1))))
             )))))

Upvotes: 0

Views: 761

Answers (3)

WorBlux
WorBlux

Reputation: 1413

Here's a general method of doing monte-carlo it accepts as arguments the number of iterations, and a thunk (procedure with no arguments) that should return #t or #f which is the experiment to be run each iteration 

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

Now it's just a mater of writing the specific experiment

You could write in your experiment where experiment is invoked in monte-carlo, but abstracting here gives you a much more flexible and comprehensible function. If you make a function do too many things at once it becomes hard to reason about and debug. 

(define (marble-experiment)
 (let ((x ...)  ;;assuming you can come up with 
       (y ...)) ;;a way to get a random x between 0 and 1
                ;;with sufficient granularity for your estimate)
  (< (sqrt (+ (* x x) (* y y))) 1)))

(define pi-estimate
  (* 4 (monte-carlo 1000 marble-experiment)))

Upvotes: 0

user448810
user448810

Reputation: 17851

I wrote a solution to this problem at my blog; the inner function is called sand because I was throwing grains of sand instead of marbles:

(define (pi n)
    (define (sand?) (< (+ (square (rand)) (square (rand))) 1))
    (do ((i 0 (+ i 1)) (p 0 (+ p (if (sand?) 1 0))))
        ((= i n) (exact->inexact (* 4 p (/ n))))))

This converges very slowly; after a hundred thousand iterations I had 3.14188. The blog entry also discusses a method for estimating pi developed by Archimedes over two hundred years before Christ that converges very quickly, with 27 iterations taking us to the bound of double-precision arithmetic.

Upvotes: 1

C. K. Young
C. K. Young

Reputation: 222973

Your parentheses are all messed up, and your argument order for < is wrong. Here's how the code should look like after it's corrected:

(define (monte-carlo-sampling n)
  (let ((x (- (* 2 (random)) 1))
        (y (- (* 2 (random)) 1)))
    (cond ((= n 0)
           0)
          (else
           (cond ((< (sqrt (+ (square x) (square y))) 1)
                  (+ 1 (monte-carlo-sampling (- n 1))))
                 (else
                  (monte-carlo-sampling (- n 1))))))))

This returns the number of hits. You'd have to convert the number of hits into a pi estimate using an outer function, such as:

(define (estimate-pi n)
  (* 4 (/ (monte-carlo-sampling n) n)))

Here's how I'd write the whole thing, if it were up to me:

(define (estimate-pi n)
  (let loop ((i 0)
             (hits 0))
    (cond ((>= i n)
           (* 4 (/ hits n)))
          ((<= (hypot (sub1 (* 2 (random)))
                      (sub1 (* 2 (random)))) 1)
           (loop (add1 i) (add1 hits)))
          (else
           (loop (add1 i) hits)))))

(Tested on Racket, using the definition of hypot I gave in my last answer. If you're not using Racket, you have to change add1 and sub1 to something appropriate.)

Upvotes: 2

Related Questions