;-------------------------------------------------------------------------------------------------
; Ex.5, Gordon et al., Prob.Prog., 2014 --- E[i,d,g,s,l] = (0.30, 0.40, ...) ---
; PCM 2014/08/02
;-------------------------------------------------------------------------------------------------                                                                                                                                                

(define no_of_samples 60000)
(define no_of_vars 5)           ; I(=Intelligence), D(=Difficulty), G(=Grade), S(=SAT), L(=Letter)
(define math_expect 0.00576)    ; mathematical expectation  (Gordon et al., 2014)

(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)
  (define i (if (flip 0.3) 1 0))                                  ; generative model
  (define d (if (flip 0.4) 1 0))                                  ; generative model
  (define g                                                       ; generative model
    (cond ((and (= (! i) 1) (= (! d) 1)) (if (flip 0.70) 1 0))
          ((and (= (! i) 1) (= d 1))     (if (flip 0.95) 1 0))
          ((and (= i 1)     (= (! d) 1)) (if (flip 0.10) 1 0))
          (else                          (if (flip 0.50) 1 0))))
  (define s                                                       ; generative model
    (if (= (! i) 1)
        (if (flip 0.05) 1 0) 
        (if (flip 0.80) 1 0))) 
  (define l                                                       ; generative model
    (if (= (! g) 1)
        (if (flip 0.10) 1 0) 
        (if (flip 0.60) 1 0)))
  (list i d g s l))                                               ; sampled value

(define (my_return)
  (let* ((time_start (get-time))
         (header 
          (display "Ex.5, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/08/02 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment2 (display "sample size = " no_of_samples))
         (comment1 
          (display "E[I, D, G, S, L] = (0.3, 0.4, ... )"))
         (samples (repeat no_of_samples take-a-sample))
         (dummy_value (hist samples "sample estimate of P(I, D, G, S, L)"))
         (comment3 (display "sample estimator of E[I, D, G, S, L] = " 
                            (expected_values samples no_of_vars)))
                  
         (comment1                            ; mathematical specification
          (display 
           "specs: P(i0, d0, g1, s1, l0) = 0.00576   (Gordon et al., 2014)"))

         (P_i1_d0_g1_s1_l0                                   ; P(I=1,D=0,G=1,S=1,L=0)
          (mean (map (lambda(sample)
                       (if (and (= (first sample) 1)         ; I(=Intelligence)
                                (= (second sample) 0)        ; D(=Difficulty)
                                (= (third sample) 1)         ; G(=Grade)
                                (= (fourth sample) 1)        ; S(=SAT)
                                (= (fifth sample) 0))        ; L(=Letter
                           1 0)) samples)))

         (comment4 (display "sample estimator of P(I=1,D=0,G=1,S=1,L=0)"  P_i1_d0_g1_s1_l0))
         (commentX(display "|deviation| = " (abs (- math_expect P_i1_d0_g1_s1_l0))))
         (time_stop (get-time)))
    (display "computation time in sec =" (/ (- time_stop time_start) 1000))))

(my_return) ; ==> E[ I, D, G, S, L] = (0.2962 0.4027333333333 0.639866666666 0.2719333333333 0.4229)

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