;------------------------------------------------------------------------------------------------- ; Ex.6b, Gordon et al., Prob.Prog., 2014 --- E[ G | L = 1] = 0.9136038186157518 --- ; ; PCM 2014/08/02 ;------------------------------------------------------------------------------------------------- (define no_of_samples 20000) (define no_of_vars 1) ; G(=Grade) (define math_expect 0.9136038186157518) ; obtained by FIGARO's variable elimination (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 (= l 1) ; observational constraint, evidence: L = 1 (Letter = weak) (list g) ; sampled value: Grade (take-a-sample))) (define (my_return) (let* ((time_start (get-time)) (header (display "Ex.6b, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/08/02 ***")) (line (display "--------------------------------------------------------------------------------")) (comment1 (display "sample size = " no_of_samples)) (samples (repeat no_of_samples take-a-sample)) (dummy_value (hist samples "sample estimate of P(G | L = 1)")) (sample-based-estimator (first (expected_values samples no_of_vars))) (comment2 (display "by FIGARO's exact variable elimination algorithm:" "parameter theta E[G | L = 1] = " math_expect)) (comment3 (display "by CHURCH's approximate rejection sampling:" "sample-based estimator E[G | L = 1] = " 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[ G | L = 1] = 0.91305 (result of one computer run with 20000 samples) ;-------------------------------------------------------------------------------------------------