;-------------------------------------------------------------------------------------------------
; 3.1, Lee&Wagenmakers, Bayesian Cognitive Modeling ---  E(theta | k=5,n=10) = 0.5
;
; PCM 2014/08/10
;-------------------------------------------------------------------------------------------------

(define no_of_samples 20000)
(define no_of_vars 1)                              ; theta
(define math_expect 0.5)                           ; mathematical expectation
(define list_of_vars '(theta))                     ; list of vars

(define (expected_value nth samples)
  (let* ((list_of_vals (map (lambda (sample) (list-ref sample nth)) samples))
         (my_mean (mean list_of_vals))
         (dummy (display "sample-based mean of " (list-ref list_of_vars nth) my_mean))
         (title_of_density (string-append "P(" (list-ref list_of_vars nth) " | k=5, n=10)"))
         (dummy (density list_of_vals title_of_density))
         (list_of_binned_vals (map (lambda(sample) (/ (round (* sample 10)) 10)) list_of_vals))
         (dummy (hist list_of_binned_vals "P(theta | k=5, n=10)"))
         (line (display "---------------------------------------------------------------------")))
    my_mean))

(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 theta (beta 1 1))                        ; generative model 
  (define k (binomial theta 10))                   ; generative model 
  (if (= k 5)                                      ; observational constraint, evidence
      (list theta)                                 ; sampled value  
      (take-a-sample)))                            ; recursive call

(define (my_return)                                
  (let* ((start_time (get-time))
         (header 
          (display "3.1 Inferring a Rate, Lee&Wagmk., BCM  *** CHURCH-code by PCM 2014/08/10 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------"))
         (comment1 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample))       
         (sample_based_expect (first(expected_values samples no_of_vars)))
         (comment2 
          (display 
           "specs: E(theta | k=5,n=10) =" math_expect " - published by Lee&Wagenmakers, 2013"))
         (comment3 (display "sample-based estimator of E(theta | k=5, n=10) ="sample_based_expect))
         (comment4 (display "|deviation| = " (abs (- math_expect sample_based_expect))))
         (computation_time (/ (- (get-time) start_time) 1000)))
    (display "computation time" computation_time "sec")))

(my_return) ; ==> 0.50           = exemplary result of one program run with n random experiments

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