;-------------------------------------------------------------------------------------------------
; Ex.4, Gordon et al., Prob.Prog., 2014 --- E[B] = 2/3 ---
; PCM 2014/07/30
;-------------------------------------------------------------------------------------------------

(define no_of_samples 10000) 
(define no_of_vars 1)                  ; B (= no of vars of conditional distribution) 
(define math_expect (/ 2 3))           ; mathematical expectation

(define (expected_value nth samples no_of_samples)
  (mean (map (lambda (sample) (list-ref sample nth)) samples)))

(define (expected_values samples vars)
  (let ((no_of_samples (length samples)))
    (if (= vars 0) '()
        (cons (expected_value (- no_of_vars vars) samples no_of_samples)
              (expected_values samples (- vars 1))))))   

(define (! b)        ; invert b: 0 -> 1        
  (abs (- b 1)))     ; invert b: 1 -> 0  

(define (take-a-sample)
  (take-a-sample-with 1))

(define (take-a-sample-with b)
  (define c (flip))                    ; generative model
  (if c 
      (take-a-sample-with (! b)) 
      (list b)))                       ; sampled value, 'while'                

(define (my_return)                                
  (let* ((header 
          (display "Ex.4, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/07/30 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1                            ; mathematical specification
          (display 
           "specs: E[B] = 1/2 + 1/2**3 + 1/2**5 + ... = 2/3"))
         (comment2 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample))       
         (dummy_value (hist samples "sample-based estimate of P(B)")) 
         (sample_based_expect (first (expected_values samples no_of_vars)))
         (comment3 (display "sample-based estimator of E[B] = " sample_based_expect)))
    (if (not (null? math_expect))
        (display "deviation from math. expectation = " (abs (- math_expect sample_based_expect)))
        "OK")))

(my_return) ; ==> 0.6666           = exemplary result of one program run with n random experiments

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