;----------------------------------------------------------------------------------------------------
; Ex.1a, Gordon et al., Prob.Prog., 2014 --- E[C1,C2)| true] = E[C1,C2] = (1/2, 1/2) ---
; PCM 2014/07/22
;----------------------------------------------------------------------------------------------------

(define no_of_samples 1000)                         ; number of samples (= random experiments)                                                                     

(define (expected_value selector list_of_lists)     ; j-th mean of n multivariate samples                                             
  (mean (map selector list_of_lists)))              ; j = 1,...,m;  m=2;  i = 1,...,n                                           

(define (take-a-sample)                             ; generation of a sample                                                                                                    
  (define c1 (if (flip) 1 0))                       ; sampling of random variable C1                                            
  (define c2 (if (flip) 1 0))                       ; sampling of random variable C2                                           
  (pair c1 c2))                                     ; returned value, sample                                                          ; condition, evidence, observation                             

(define (my_return)                                 ; returns results of the set of experiments                              
  (let* ((header 
          (display "Ex.1a, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/07/22 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1                             ; mathematical specification
          (display 
           "specification: E[C1,C2] = 1/4*((0,0)+(0,1)+(1,0)+(1,1)) = 1/4*(2,2) = (1/2,1/2))"))
         (comment2 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample))     ; 'samples' is bound to list of samples
         (dummy_value (hist samples "sample-based estimate of P(C1, C2)"))   ; histogram of samples
         (comment3 (display "sample-based estimator of E[C1, C2] = ")))
    (pair(expected_value first samples)     ; E[C1 | true] = E[C1] = 1/2     ; first component
         (expected_value rest  samples))))  ; E[C2 | true] = E[C2] = 1/2     ; second component 

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

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