;-------------------------------------------------------------------------------------------------
; Ex.6d, Bayesian Network 'Student Model' ---  E[I| D=0, L=1, S=0]= 0.029197080291970802  ----
;        with D=0 (difficult), L=1 (weak letter), S=0 (low SAT score)
;
; PCM 2014/08/06
;-------------------------------------------------------------------------------------------------

(define no_of_samples 20000)
(define no_of_vars 1)                    ; i (= intelligence)
(define math_expect 0.029197080291970802); by FIGARO's exact variable elimination inference algor.

(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)))
  (if  (and (= d 0)(= l 1)(= s 0))     ; condition, evidence, observation
       (list i)                        ; sampled value
       (take-a-sample)))               ; resampling if failure of observational condition

(define (my_return)
  (let* ((time_start (get-time))
         (header 
          (display "Ex.6d, BN 'Student Model'   *** CHURCH-code by PCM 2014/08/06 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample))
         (dummy_value (hist samples "sample-based estimate of P(I | D = 0, L = 1, S = 0)"))
         (sample-based-estimator (first (expected_values samples no_of_vars)))
         (comment2 (display "by FIGARO's exact variable elimination algorithm:"
                            "parameter theta E[I | D = 0, L = 1, S = 0] = " math_expect))
         (comment3 (display "by CHURCH's approximate rejection sampling:"
                            "sample-based estimator E[I | D = 0, L = 1, S = 0] = " 
                            sample-based-estimator))
         (comment4(display "|deviation| = " (abs (- math_expect sample-based-estimator))))
         (time_stop (get-time)))
    (display "computation time in sec =" (/ (- time_stop time_start) 1000))))
  

(my_return) ; E[I|D=0, L=1, S=0] = 0.029197080291970802

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