;----------------------------------------------------------------------------------------------------
; Ex.7, Gordon et al., Prob.Prog., 2014 ---  E[X] = 81/6 = 13.5 ---
; PCM 2014/07/25
;----------------------------------------------------------------------------------------------------

(define no_of_samples 30000)  
(define no_of_vars 1)                              ; X
(define math_expect (/ 81 6))                      ; 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)
  (take-a-sample-with 0))                                                  ; generative model

(define (take-a-sample-with x)                     
  (define coin (flip 0.5))                                                 ; generative model
  (cond ((= x 0)(if coin (take-a-sample-with  1) (take-a-sample-with  2))) ; generative model
        ((= x 2)(if coin (take-a-sample-with  5) (take-a-sample-with  6))) ; generative model
        ((= x 1)(if coin (take-a-sample-with  3) (take-a-sample-with  4))) ; generative model
        ((= x 3)(if coin (take-a-sample-with  1) (take-a-sample-with 11))) ; generative model
        ((= x 4)(if coin (take-a-sample-with 12) (take-a-sample-with 13))) ; generative model
        ((= x 5)(if coin (take-a-sample-with 14) (take-a-sample-with 15))) ; generative model
        ((= x 6)(if coin (take-a-sample-with 16) (take-a-sample-with  2))) ; generative model
        (else (list x))))                                                  ; cond. sampled value

(define (my_return)                                
  (let* ((header 
          (display "Ex.7, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/07/25 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1                                                  ; mathematical specification
          (display 
           "specs: E[X] = 1/6*(11 + 12 + 13 + 14 + 15 + 16) = 81/6 = 13.5"))
         (comment2 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample))       
         (dummy_value (hist samples "sample-based estimate of P(X)")) 
         (sample_based_expect (first (expected_values samples no_of_vars))) ; E[X] = 13.5 
         (comment3 (display "sample-based estimator of E[X] = " sample_based_expect))) 
    (if (not (null? math_expect))
        (display "deviation from math. expectation = " (abs (- math_expect sample_based_expect)))
        "OK")))

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

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