;-------------------------------------------------------------------------------------------------
; Ex.8, Gordon et al., Prob.Prog., 2014 ---  Bayesian Skill Rating (TrueSkill model) ---
; PCM 2014/08/10
;-------------------------------------------------------------------------------------------------

(define no_of_samples 20000)                                        ; number of random experiments
(define no_of_vars 5)  
(define list_of_vars '(skillA skillB skillC skill_precision perf_precision))
(define math_expect '(105.7 100.0 94.3 0.11 0.11))              ; published by Gordon et al., 2014

(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 "mean of " (list-ref list_of_vars nth) my_mean))
         (title_of_density (string-append "P(" (list-ref list_of_vars nth) " | Outcomes)"))
         (dummy (density list_of_vals title_of_density))
         (line (display "---------------------------------------------------------------------")))
    my_mean))

(define (expected_values samples vars)
  (if (= vars 0) '()
      (cons (expected_value (- no_of_vars vars) samples)
            (expected_values samples (- vars 1)))))

(define (take-a-sample)
  ; standard deviations
  (define shape 8.0)              ; "shape"-(alpha)-Parameter of gamma(shape, scale)
  (define scale 1.0)              ; "scale"-(beta)-Parameter  of gamma(shape, scale)
  (define skill_sd (gamma shape scale)) ; mean = shape*scale = 8 ; std = sqrt(shape)*scale = 2.83
  (define perf_sd  (gamma shape scale)) ; mean = shape*scale = 8 ; std = sqrt(shape)*scale = 2.83
  ; skill model
  (define skillA (gaussian 100.0  skill_sd))  ; generative model
  (define skillB (gaussian 100.0  skill_sd))  ; generative model
  (define skillC (gaussian 100.0  skill_sd))  ; generative model
  ; first game: A vs B, A won
  (define perfA1 (gaussian skillA perf_sd))   ; generative model
  (define perfB1 (gaussian skillB perf_sd))   ; generative model
  ; second game: B vs C, B won
  (define perfB2 (gaussian skillB perf_sd))   ; generative model
  (define perfC2 (gaussian skillC perf_sd))   ; generative model
  ; third game: A vs C, A won
  (define perfA3 (gaussian skillA perf_sd))   ; generative model
  (define perfC3 (gaussian skillC perf_sd))   ; generative model
  ; rejection sampling
  (if (and (> perfA1 perfB1)(> perfB2 perfC2)(> perfA3 perfC3)); observations (= game performances)
      (list skillA skillB skillC (/ 1.0 skill_sd)(/ 1.0 perf_sd)); samples (skills, precisions)
      (take-a-sample)))

(define (my_return)
  (let* ((start_time (get-time))
         (header 
          (display "Ex.8, Gordon et al., Prob.Prog., 2014 *** CHURCH-code by PCM 2014/08/10 ***"))
         (line (display "---------------------------------------------------------------------"))
         (comment1 (display "specs: " math_expect " - published by Gordon et al., 2014"))
         (comment2 (display "sample size = " no_of_samples))
         (line (display "---------------------------------------------------------------------"))
         (samples (repeat no_of_samples take-a-sample))
         (sample_based_means(expected_values samples no_of_vars))
         (comment3 (display "mathematical expectation =" math_expect))
         (comment4 (display "sample-based means = " sample_based_means))
         (sum_of_abs_dev 
          (sum (map (lambda (x y) (abs(- x y))) math_expect sample_based_means)))
         (dummy3 (display "Sum of |deviation| =" sum_of_abs_dev )) ; ==> 0.01384601410834941
         (computation_time (/ (- (get-time) start_time) 1000)))
    (display "computation time" computation_time "sec")))
    
;result of exemplary run: E[skillA skillB skillC skill_precision perf_precision | outcomes] 
; = (105.69271200421967 99.90985708384883 94.1775550319377 0.13968936710191396 0.14589155872673595)
(my_return)
;-------------------------------------------------------------------------------------------------