;-------------------------------------------------------------------------------------------------
; 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

;-------------------------------------------------------------------------------------------------