;------------------------------------------------------------------------------------------------- ; file: 'My_Pi' --- P(x^2+y^2 <= 1 | 0<=x<=1,0<=x<=1) = pi/4 --- ; PCM 2014/08/10 ;------------------------------------------------------------------------------------------------- (define no_of_samples 2000000) ; first hundred digits of pi (define pi 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679) (define math_expect (/ pi 4.0)) ; P(x^2+y^2 <= 1 | 0<=x<=1,0<=x<=1) = pi/4 (define epsilon 0.001) ; epsilon for Hoeffding bound (Koller & Friedman, A.2.2) (define (hoeffding_bound eps m) ; Koller & Friedman, A.2.2 (exp (* -2.0 m (* eps eps)))) ; p = Bernoulli parameter, m = no_of_samples, eps = epsilon (define (take-a-sample) (define x (uniform 0.0 1.0)) ; generative model (one side of unit square) (define y (uniform 0.0 1.0)) ; generative model (one side of unit square) (define (inside_unit_circle? x y) ; indicator function for 'inside' unit circle (if (<= (+ (* x x)(* y y)) 1.0) 1 ; sampled value, if 'inside' unit-circle 0)) ; sampled value, if 'outside' unit-circle (inside_unit_circle? x y)) (define (my_return) (let* ((t_start (get-time)) (header (display "file: 'My_Pi' *** CHURCH-code by PCM 2014/08/10 ***")) (line (display "--------------------------------------------------------------------------------")) (comment1 ; mathematical specification (display "math. expectation p: P(x^2 + y^2 <= 1 | 0<=x<=1, 0<=y<=1) = pi/4 =" math_expect)) (comment2 (display "sample size = " no_of_samples)) (samples (repeat no_of_samples take-a-sample)) (dummy_value (hist samples "P(x^2+y^2 <= 1 | 0<=x<=1, 0<=x<=1)")) (#inside (sum samples)) ; number of particles inside unit circle (#total (length samples)) ; number of particles inside unit square (sample_based_expect (/ #inside #total)) (comment3 (display "sample-based estimator TD (proportion of successes) of " "P((x,y) is inside unit circle | (x,y) is inside unit square) = (#inside / #total) = " sample_based_expect)) (deviation (- math_expect sample_based_expect)) (comment4 (display "deviation = (math. expectation - sample_based_expect) = " deviation)) (comment5 (display "pi = " pi)) (sample_based_pi (* 4.0 sample_based_expect)) (comment6 (display "sample-based estimator of pi = " sample_based_pi)) (comment7 (display "epsilon for Hoeffding bound = " epsilon)) (hoeffding (hoeffding_bound epsilon no_of_samples)) (comment8 (if (>= deviation 0) (display "inequality for positive excess:" " P(TD > p + eps) <= exp(-2*M*eps^2)) = " hoeffding "P(#inside/#total > pi/4 + eps) <= exp(-2*M*eps^2) = " hoeffding) (display "inequality for negative excess:" " P(TD < p - eps) <= exp(-2*M*eps^2)) = " hoeffding "P(#inside/#total < pi/4 - eps) <= exp(-2*M*eps^2) = " hoeffding))) (t_stop (get-time))) (display "computation time = " (/ (- t_stop t_start) 1000) "sec"))) (my_return) ; ==> 0.785264 = result of one program run with 2.000.000 random trials ;-------------------------------------------------------------------------------------------------