;--------------------------------------------------------------------------------------------------
; Ex.1b, Gordon et al., Prob.Prog., 2014 --- E[C1,C2)| c1 v c2] = (2/3, 2/3) ---
; PCM 2014/11/15
;--------------------------------------------------------------------------------------------------

(define (hoeffding_invers eps prob_bound)           ; n >= (ln(2) - ln(bound))/(2*eps^2)
  (round (/ (- (log 2) (log prob_bound)) (* 2.0 eps eps))))

(define no_of_samples (hoeffding_invers 0.01 0.05)) ; 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                                           
  (if (or (= c1 1) (= c2 1))                        ; condition, evidence, observation   
      (pair c1 c2)                                  ; returned value, sample                             
      (take-a-sample)))                             ; rejection sampling

(define (my_return)                                 ; returns results of the set of experiments
  (let* ((start_time (get-time))
         (header 
          (display "Ex.1b, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/11/15 ***"))
         (line
          (display 
           "------------------------------------------------------------------------------------"))
         (comment1                             ; mathematical specification
          (display 
           "specs: E[C1,C2|C1=1 v C2=1] = 0*(0,0)+1/3((0,1)+(1,0)+(1,1)) = 1/3*(2,2) = (2/3,2/3)"))
         (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 | C1=1 v C2=1)"))
         (comment3 (display "sample-based estimator of E[(C1, C2) | C1=1 v C2=1] = "))
         (comment4 (display 
                    (pair (expected_value first samples)      ; E[C1 | c1 v c2 ] = 2/3 -- first  component --
                          (expected_value rest  samples))))   ; E[C2 | c1 v c2 ] = 2/3 -- second component --
         (computation_time (/ (- (get-time) start_time) 1000)))
    (display "computation time in sec " computation_time)))
    
(my_return) ; ==> (0.6697 . 0.6665) = exemplary result of one program run with n random experiments

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