;----------------------------------------------------------------------------------------------------------
; 4.1, Lee&Wagenmakers, Bayesian Cognitive Modeling    ---   Infering a mean and a standard deviation   ---
;                                                      ---            Importance Sampling               ---
; PCM 2014/10/12
;----------------------------------------------------------------------------------------------------------
(define (stddev lst_x)
  (define (square x) (* x x))
  (sqrt (- (/ (apply + (map square lst_x)) (length lst_x)) (square (mean lst_x)))))

(define (expected_values samples weights vars)
  (let ((no_of_samples (length samples)))
    (if (= vars 0) '()
        (cons (expected_value (- no_of_post_vars vars) samples weights)
              (expected_values samples weights (- vars 1))))))

(define (expected_value nth samples weights)
  (let* ((list_of_vals 
          (map (lambda (sample) (list-ref sample nth)) samples))
         (list_of_weighted_vals 
          (map (lambda (sample weight) 
                 (* (list-ref sample nth) weight)) samples weights))
         (my_mean (sum list_of_weighted_vals))
         (dummy (display "sample-based mean of " 
                         (list-ref list_of_post_vars nth) my_mean "------------------"))
         (my_stddev (stddev list_of_weighted_vals))
         (dummy (display "sample-based stddev of " 
                         (list-ref list_of_post_vars nth) my_stddev "------------------"))
         (title_of_density 
          (string-append "P(" (list-ref list_of_post_vars nth) " | xs = 1.1, 1.9, 2.3, 1.8) " ))
         (dummy (density list_of_weighted_vals title_of_density))
         (list_of_binned_vals (map (lambda(sample) (/ (round (* sample 10)) 10)) list_of_weighted_vals))
         (dummy (hist list_of_binned_vals title_of_density))
         (line (display "------------------------------------------------------------------------------")))
    my_mean))

(define (log-likelihood mu sigma lst_x)  ; Fahrmeier et al, Statistik, 4/e, 2003, p.377
  (define (square x) (* x x))
  (define n (length lst_x))
  (define my_pi 3.141592654)
; (define factor (* n (log (sqrt (* 2 my_pi)))))   ; - n*ln(sqrt(2pi)) ; cancels out in the ratio
  (define ssq (apply + (map (lambda (x) (square (- x mu))) lst_x)))
  (- 0
;    factor                                        ; - n*ln(sqrt(2pi)) ; cancels out in the ratio
     (* n (log sigma))                             ; - n*ln(sigma)
     (* (/ 1 (* 2 (square sigma)))
        ssq)))
;----------------------------------------------------------------------------------------------------------

(define no_of_samples 70000)
(define math_expect 1.774)                                    ; expectation obtained from OpenBUGS
(define list_of_post_vars '(mu sigma lambda accept_ratio))    ; list of aposterior vars
(define no_of_post_vars (length list_of_post_vars))           ; mu, sigma, lambda
(define xs_data '(1.1 1.9 2.3 1.8))                           ; data: 1.1, 1.9, 2.3, 1.8
(define no_of_data (length xs_data))                          ; n 
(define mean_x (mean xs_data))
(define stddev_x (stddev xs_data))
(define log-max-likelihood (log-likelihood mean_x stddev_x xs_data))

(define (take-a-sample)  
  (define mu (gaussian 0 31.6227766))        ; prior of mu, precision = .001 => s = 31.62...
  (define sigma (uniform 0 10))              ; prior of sigma
  (define my_lambda (/ 1 (expt sigma 2)))
  (define accept_ratio
    (exp (- (log-likelihood mu sigma xs_data) log-max-likelihood)))
  (list mu sigma my_lambda accept_ratio))    ; weighted sample values 

(define (my_return)                                
  (let* ((start_time (get-time))
         (header 
          (display 
           "4.1 - 'Inferring a mean and std.dev.' Lee&Wagmk., BCM  *** CHURCH-code by PCM 2014/10/12 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------------------"))
         (comment1 (display "sample size = " no_of_samples))
         (samples (repeat no_of_samples take-a-sample)) 
         (sum_of_weights (apply + (map (lambda (sample) (last sample)) samples)))
         (dummy (display "sum_of_weights " sum_of_weights))
         (weights (map (lambda (sample) (/ (last sample) sum_of_weights)) samples))
         (sum_of_normalized_weights (sum weights))
         (dummy (display "sum_of_normalized_weights " sum_of_normalized_weights))
         (sample_based_expect (first (expected_values samples weights no_of_post_vars)))
         (comment2 
          (display 
           "specs: E(mu, sigma | xs = 1.1, 1.9, 2.3, 1.8) =" 
           math_expect " - published by Lee&Wagenmakers, 2013"))
         (comment3 (display "sample-based estimator of E(mu, sigma | xs = 1.1, 1.9, 2.3, 1.8) ="
                            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) ;            = exemplary result of one program run with n random experiments

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