#| 

  AR(1)  Racket Gamble.

  From Stan Users Guide, section 3.1 (page 51) "Autoregressive Models"
  https://mc-stan.org/docs/2_19/stan-users-guide-2_19.pdf
  (page 48f)

  Mathematica:
  """
  eproc = EstimatedProcess[temp, ARProcess[1]]
  -> ARProcess[0.731788, {0.157416}, 0.118305]
  """

var : alpha
1.1090322845052043: 0.708590044921784
0.5491458279911212: 0.29140995507821604
0.685890970444901: 2.0593060111469593e-20
0.44631079038563315: 1.5713639007088746e-25
...
1.6448808668533812: 3.6164913541496907e-300
3.1678604032510798: 2.6138188566346616e-300
4.206989733279477: 3.1733067212797045e-302
-4.126772034141277: 0.0
mean: 0.9458757973635339

var : beta
-0.37045143739471253: 0.708590044921784
0.4046769344729959: 0.29140995507821604
0.4508831857811233: 2.0593060111469593e-20
0.24893528291429767: 1.5713639007088746e-25
...
0.9194451359428357: 3.6164913541496907e-300
0.6749488766280096: 2.6138188566346616e-300
0.7468932705094318: 3.1733067212797045e-302
0.3124909509366822: 0.0
mean: -0.1445713133688929

var : sigma
0.3675208514860001: 0.708590044921784
0.4019351517816607: 0.29140995507821604
0.5988646460886968: 2.0593060111469593e-20
0.9533285559572817: 1.5713639007088746e-25
...
12.871286474814879: 3.6164913541496907e-300
2.6607404510467383: 2.6138188566346616e-300
5.383797322375248: 3.1733067212797045e-302
1.6870507302202808: 0.0
mean: 0.3775495211892068



  This is a port of my WebPPL model ar1.wppl 
 
  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 ys '(
             0.705429 1.43062 0.618161 0.315107 1.09993 1.49022 0.690016 0.587519
             0.882984 1.0278 0.998615 0.878366 1.17405  0.532718 0.486417
             1.13685 1.32453 1.3661 0.914368 1.07217 1.1929 0.418664
             0.889512 1.47218 1.13471 0.410168 0.639765 0.664874 1.12101
             1.22703 -0.0931769 0.4275 0.901126 1.01896 1.27746 1.17844
             0.554775 0.325423 0.494777 1.05813 1.10177 1.11225 1.34575
             0.527594 0.323462 0.435063 0.739342 1.05661 1.42723 0.810924
             0.0114801 0.698537 1.13063 1.5286 0.968813 0.360574 0.959312
             1.2296 0.994434 0.59919 0.565326 0.855878 0.892557 0.831705
             1.31114 1.26013 0.448281 0.807847 0.746235 1.19471 1.23253
             0.724155 1.1464 0.969122 0.431289 1.03716 0.798294 0.94466
             1.29938 1.03269 0.273438 0.589431 1.2741 1.21863 0.845632
             0.880577 1.26184 0.57157 0.684231 0.854955 0.664501 0.968114
             0.472076 0.532901 1.4686 1.0264 0.27994 0.592303 0.828514
             0.625841
             ))

(define (model)

  (; enumerate
   ; rejection-sampler
   importance-sampler
   ; mh-sampler ; #:transition (slice) 
   
   (define alpha (normal 2 2))
   (define beta  (normal 2 2))
   ; (define sigma (gamma 2 2))
   (define sigma (sample (clip-distx (gamma-dist 2 2) 0 10000)))
   
   ; (observe/fail (> sigma 0))
   
   (define (y i)
     (if (= i 0)
         (normal-dist (+ alpha beta) sigma)
         (normal-dist (+ alpha (* beta (sample (y (sub1 i))))) sigma)))
   
   (for ([i (length ys)])
     (observe-sample (y i) (list-ref ys i)))
   
   (list alpha
         beta
         sigma
         )

   )
  )

(show-marginals (model)
                (list "alpha"
                      "beta"
                      "sigma"
                        )
                  #:num-samples 1000
                  #:truncate-output 4
                  ; #:skip-marginals? #t
                  ; #:credible-interval 0.94
                  ; #:show-stats? #t
                  ; #:show-histogram? #t
                  )