;------------------------------------------------------------------------------------------------- ; Ex.3, Gordon et al., Prob.Prog., 2014 --- E[Count| c1 v c2] = 4/3 ; PCM 2014/07/28 ;------------------------------------------------------------------------------------------------- (define no_of_samples 10000) (define no_of_vars 1) ; for Count (define math_expect (/ 4 3)) (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 (take-a-sample) (define c1 (if (flip) 1 0)) ; generative model (define c2 (if (flip) 1 0)) ; generative model (define count (+ c1 c2)) ; generative model (if (not (or (= c1 1)(= c2 1))) ; negated observational constraint (take-a-sample) ; the recursive call substitutes the 'while' of PROB (list count))) ; return value of 'take-a-sample' (define (my_return) (let* ((header (display "Ex.3, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/07/28 ***")) (line (display "--------------------------------------------------------------------------------")) (comment1 ; mathematical specification (display "specs: E[Count|c1 v c2] = 0*0+1/3*(1+1+2) = 1/3(4) = 4/3")) (comment2 (display "sample size = " no_of_samples)) (samples (repeat no_of_samples take-a-sample)) (dummy_value (hist samples "sample-based estimate of P(Count | c1 v c2)")) (sample_based_expect (first (expected_values samples no_of_vars))) (comment3 (display "sample-based estimator of E[Count| c1 v c2] = " sample_based_expect))) (if (not (null? math_expect)) (display "deviation from math. expectation = " (abs (- math_expect sample_based_expect))) "OK"))) (my_return) ; ==> 1.333 = exemplary result of one program run with n random experiments ;-------------------------------------------------------------------------------------------------