#| 

  Bayesian linear regression in Racket Gamble.

  From BLOG example/bayesian-linear-regression.blog
  """
  A model for Bayesian linear regression
 
  @author yiwu
  @since 2015-07-02
  """

  This is a port of my WebPPL model bayesian_linear_regression.wppl

  Output:

var : y4
9.10161212766774: 0.12086945142623648
9.309931784525654: 0.11172337584984503
10.91709573383832: 0.08399858458426669
8.538020373246502: 0.07350715263883331
10.446309707341985: 0.058466516795837316
...
-4.595764056143208: 3.8947e-320
-0.35748359804252006: 1.891e-320
-5.445626426157869: 1.374e-321
1.4976566327198673: 1.6e-322
-31.7402435394563: 0.0
mean: 9.517461577263164
Min: 1.875802431135524 Mean: 9.890913853575805 Max: 15.869430346719128 Variance: 2.6406858871075447 Stddev: 1.6250187343866362
Credible interval (0.84): 7.677174579925393..12.039605736512556

var : m
1.8088122207250599: 0.12086945142623648
1.794759193331004: 0.11172337584984503
2.011225398087395: 0.08399858458426669
1.7005492805632993: 0.07350715263883331
1.8941054593952003: 0.058466516795837316
...
-0.9825895084680619: 3.8947e-320
-0.6124958647725994: 1.891e-320
-1.5331560219667937: 1.374e-321
0.38757089020558605: 1.6e-322
-6.4866772864873825: 0.0
mean: 1.832840272267902
Min: 0.11168098987964598 Mean: 2.0077984987338753 Max: 3.33280868249568 Variance: 0.11898522310637583 Stddev: 0.3449423475109657
Credible interval (0.84): 1.5674895399954556..2.5043307338932594

var : b
-0.15167427216153415: 0.12086945142623648
0.25119431751820104: 0.11172337584984503
0.14777102051014587: 0.08399858458426669
0.3606963570368861: 0.07350715263883331
0.32746266170015453: 0.058466516795837316
...
0.5451695804274529: 3.8947e-320
2.575498779107246: 1.891e-320
2.683243328208245: 1.374e-321
-0.41522319556717874: 1.6e-322
1.2085755471317516: 0.0
mean: 0.026415834611575255
Min: -2.5654844426053542 Mean: -0.27012904887281086 Max: 2.3807112665192767 Variance: 0.38203349811592596 Stddev: 0.6180885843598197
Credible interval (0.84): -0.8181855111657625..0.8108726072452849

var : sigma2
0.534687138917272: 0.12086945142623648
0.6076264067937871: 0.11172337584984503
0.5022786601326074: 0.08399858458426669
0.5765133248441577: 0.07350715263883331
0.5035790951883868: 0.058466516795837316
...
0.2632500621047391: 3.8947e-320
0.16889451731359764: 1.891e-320
0.24502558965038376: 1.374e-321
0.1679517801338904: 1.6e-322
0.7181242806122738: 0.0
mean: 0.9005723899748515
Min: 0.35587473505676487 Mean: 0.7333230756929937 Max: 2.9578508551119427 Variance: 0.1717278879735853 Stddev: 0.4144006370332764
Credible interval (0.84): 0.35587473505676487..1.1667806923215733


  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Racket page: http://www.hakank.org/racket/

|#

#lang gamble

; (require gamble/viz)
(require racket)
(require "gamble_utils.rkt")


(define (bayesian-linear-regression)
  
  (; enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler

   (define m (normal 0 2))
   (define b (normal 0 2))  
   (define sigma2 (gamma 1 1))
   
   (define (y x) (normal-dist (+ (* m x) b) sigma2))

   (define ys '(0 1 4 6))
   (for ([i (range (length ys))])
     (observe-sample (y i) (list-ref ys i)))

   (define y4 (sample (y 5)))
  
   (list y4 m b sigma2)
   )
  )

(show-marginals (bayesian-linear-regression) '("y4" "m" "b" "sigma2")
                #:truncate-output 5
                #:show-stats? #t
                #:credible-interval 0.84)