;----------------------------------------------------------------------------------------------------------
; 4.1, Lee&Wagenmakers, Bayesian Cognitive Modeling --- Inference of a mean and a standard deviation  ---
;
; Church code by PCM 2016/02/03 with (factor log(likelihood)+log(priors)) and burnin-period excluded
;
;                comments are welcome: claus.moebus[at]uni-oldenburg.de
;----------------------------------------------------------------------------------------------------------

;----------------------------------------------------------------------------------------------------------
;  Bravais-Pearson correlation coefficient
;----------------------------------------------------------------------------------------------------------

(define (corr xs ys) ; xs, ys are lists of data to be correlated 
  (let* ((n (length xs))
         (ssq (lambda (x y)
                (let* ((xm (mean x))
                       (ym (mean y))
                       (sxy (sum (map (lambda(xi yi) (* xi yi)) x y))))
                  (- sxy (* n xm ym))))))
    (/ (ssq xs ys) (sqrt (* (ssq xs xs) (ssq ys ys)))))) ; Fahrmeier et al., Statistik, 2003, p.137

;----------------------------------------------------------------------------------------------------------
;  ACF(lag) (autocorrelation function)
;----------------------------------------------------------------------------------------------------------

(define (autocorr xs)
  (let* ((n (length xs))
         (lags (range 1 (round (* n 0.01)))))             ; max(lag) = 1% of n
    (list (reverse lags) 
          (reverse (map (lambda (lag) (corr (take xs (- n lag)) (drop xs lag))) lags)))))
 
;----------------------------------------------------------------------------------------------------------
;  standard deviation
;----------------------------------------------------------------------------------------------------------

(define (stddev lst_x)
  (define (square x) (* x x))
  (sqrt (- (/ (sum (map square lst_x)) (length lst_x)) (square (mean lst_x)))))

;----------------------------------------------------------------------------------------------------------
; summary statistics for posterior distributions
;----------------------------------------------------------------------------------------------------------

(define (expected_value nth samples)
  (let* ((list_of_vals (map (lambda (sample) (list-ref sample nth)) samples))
         (my_mean (mean list_of_vals))
         (dummy1 
          (display "sample-based Bayesian posterior mean of " (list-ref list_of_vars nth) my_mean))
         (my_stddev (stddev list_of_vals))
         (dummy2 
          (display "sample-based Bayesian posterior standard deviation of " 
                   (list-ref list_of_vars nth) my_stddev))
         (title_of_density 
          (string-append "P(" (list-ref list_of_vars nth) " | xs = 1.1, 1.9, 2.3, 1.8) " ))
         (acf-graph (barplot (autocorr list_of_vals)
                             (string-append "ACF(Lag) for " (list-ref list_of_vars nth))))
         (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 title_of_density))
         (line (display "---------------------------------------------------------------------")))
    my_mean))

(define (expected_values samples vars)j
  (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)))))) 

;----------------------------------------------------------------------------------------------------------
; log-likelihood and log-priors
;----------------------------------------------------------------------------------------------------------

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

(define (log-prior-mu mu sigma lst_mu) 
  (log-likelihood mu sigma lst_mu))                   ; density of N(x, mu, sigma)

(define (log-prior-sigma a b) (/ 1 (- b a)))          ; density of uniform(a, b)

;----------------------------------------------------------------------------------------------------------
; Initials
;----------------------------------------------------------------------------------------------------------

(define no_of_samples 10000)
(define thin-factor 200)                              ; must be a natural number >= 1
(define length_of_burnin (* no_of_samples .20))       ; length of burning-in period (20%)
(define empirical_mean 1.775)                         ; (1.1 + 1.9 + 2.3 + 1.8)/4 = 1.775
(define list_of_vars '(Mu Sigma Lambda log_pi_theta_e))  ; list of posterior vars (eg:...log(pi(theta|e)))
(define no_of_vars (length list_of_vars))             ; mu, sigma, lambda, P_acceptance 
(define xs_data '(1.1 1.9 2.3 1.8))                   ; data: xs = (1.1, 1.9, 2.3, 1.8)
(define mean_x (mean xs_data))
(define stddev_x (stddev xs_data))

;----------------------------------------------------------------------------------------------------------
; Probabilistic Model
;----------------------------------------------------------------------------------------------------------

(define (take-samples)  
    (mh-query 
     no_of_samples 
     thin-factor                                      ; 'lag' or 'thin-factor'
     (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)))          ; prior of precision lambda = 1/sigma^2
     (define log_pi_theta_e                           ; log of posterior model density log(pi(theta|e))
       (+ (log-likelihood mu sigma xs_data)     
          (log-prior-mu 0 31.6227766 (list mu))
          (log-prior-sigma 0 10)))
     (list mu sigma my_lambda log_pi_theta_e)         ; sampled (posterior) values  
     (factor log_pi_theta_e)))                   

;----------------------------------------------------------------------------------------------------------
; comments
;----------------------------------------------------------------------------------------------------------

(define (my_run)                                
  (let* ((start_time (get-time))
         (header 
          (display 
           "4.1 - 'Inferring a mean & stddev' Lee & Wagmk.,BCM  *** CHURCH-code by PCM 2016/02/03 ***"))
         (line
          (display 
           "--------------------------------------------------------------------------------------------"))
         (comment1 (display "sample size (after thinning, burn-in period included) = " no_of_samples))
         (comment2 (display "thin factor = " thin-factor))
         (comment3 (display "length(burn-in period) = " length_of_burnin))
         (comment4 
          (display "effective sample size (burn-in period excluded) = " 
                   (- no_of_samples length_of_burnin)))
         (prior_mu (display (density (repeat no_of_samples (lambda () (gaussian 0 31.6227766))) "P(Mu)")))
         (prior_sigma (display (density (repeat no_of_samples (lambda () (uniform 0 10))) "P(Sigma)")))
         (samples (drop (take-samples) length_of_burnin)) ; exclusion of burnin-period    
         (sample_based_expect (first(expected_values samples no_of_vars)))
         (comment5 
          (display "specs: empiricial mean = " empirical_mean))
         (comment6 (display "Bayesian estimator E(mu, sigma | xs = 1.1, 1.9, 2.3, 1.8) ="
                            sample_based_expect))
         (comment7 (display "|deviation| = " (abs (- empirical_mean sample_based_expect))))
         (computation_time (/ (- (get-time) start_time) 1000)))
    (display "computation time" computation_time "sec")))

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

(my_run) ;            = exemplary result of one program run with n random experiments

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