;-------------------------------------------------------------------------------------------------
; Ex.2, Gordon et al., Prob.Prog., 2014 ---  E[Count| c1 v c2] = 4/3 
; PCM 2014/07/27
;-------------------------------------------------------------------------------------------------

(define no_of_samples 10000)
(define no_of_vars 1)                              ; Count
(define math_expect (/ 4 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 (take-a-sample)                                                                                                                                                           
  (define c1 (if (flip) 1 0))                      ; generative model                                          
  (define c2 (if (flip) 1 0))                      ; generative model                                          
  (define count (+ c1 c2))                         ; generative model
  (if (or (= c1 1)(= c2 1))                        ; observational constraint, evidence
      (list count)                                 ; sampled value  
      (take-a-sample)))                            ; recursive call

(define (my_return)                                
  (let* ((header 
          (display "Ex.2, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/07/27 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1                            ; mathematical specification
          (display 
           "specs: E[Count|c1 v c2] = 0*0+1/3*(1+1+2) = 1/3*(4) = 4/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(Count | c1 v c2)")) 
         (sample_based_expect (first (expected_values samples no_of_vars)))
         (comment3 (display "sample-based estimator of E[Count| c1 v c2] = " sample_based_expect)))
    (if (not (null? math_expect))
        (display "deviation from math. expectation = " (abs (- math_expect sample_based_expect)))
        "OK")))

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

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