#| 

  Changepoint detection in Racket Gamble.

  From https://nishanthu.github.io/articles/ChangePointAnalysis.html

  This is a port of my WebPPL model changepoint.wppl

(13 : 0.5947518313314275)
(12 : 0.4050817841414437)
(8 : 7.712347581108273e-5)
(14 : 7.473693451959234e-5)
(11 : 1.3655718538874703e-5)
(10 : 7.318695731032216e-7)
(9 : 8.798372576300571e-8)
(15 : 3.939228410524439e-8)
(7 : 7.86960539434155e-9)
(6 : 7.44427253370055e-10)
(17 : 2.673542482135201e-10)
(16 : 1.3314698655190967e-10)
(5 : 7.02367845827889e-11)
(4 : 3.926361008211297e-11)
(18 : 2.454668378808425e-11)
(19 : 3.768068787278835e-12)
(21 : 2.133507409394572e-13)
(22 : 3.5780005907917906e-14)
(2 : 2.627262685891652e-14)
(23 : 1.891276002646139e-14)
(3 : 1.736426713279489e-14)
(24 : 8.078407507787324e-15)
(26 : 2.928997698333554e-15)
(0 : 1.6288743766039666e-15)
(1 : 1.606675039583082e-15)
(20 : 8.902028332342473e-16)
(27 : 2.338162790909014e-16)
(32 : 1.491981665204432e-16)
(31 : 1.37573445590999e-16)
(29 : 1.880206271517716e-17)
(30 : 6.5684249321707185e-18)
(25 : 6.029046843341405e-18)
(28 : 1.4783039134430062e-18)
(mean: 12.594577503490054)
Min: 12 Mean: 12.995 Max: 13 Variance: 0.004975 Stddev: 0.07053367989832943
Credible interval (0.94): 13..13



  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")

;This was generated by R:
; > c(rnorm(13, mean=0, sd=0.5), rnorm(20, mean=4, sd=0.5))
; I.e. the changepoint is at 13
; 
(define *ys* '(0.90685440 0.38219559 0.13508408 -0.96054590 0.78553871 0.06791336 -0.09750806
              -0.85431057 0.33478539 0.76388689 0.52911522 -0.41226421 1.53302740 
              5.01118838 4.75164473 3.79149212 3.97544567 4.33985561 4.52511631 4.21723600 
              4.22733898 3.80925885 3.65504091 3.25267865 3.85339052 3.40136529 4.21462117 
              3.94033196 4.00900359 4.44148942 4.23449516 4.52032000 4.41335270))

(define (model ys)

  (importance-sampler
   
   (define n (length ys))
   (define mu (avg ys))
   
   (define mu1 (uniform 0 mu))
   (define mu2 (uniform 0 mu))  
   
   (define cp (random-integer n))
   
   (define(mu0 i) (if (< i cp) mu1 mu2))
   
   (define y (for/list ([i (range n)])
               (normal-dist (mu0 i) 1))
     )
   
   ;; observe the values in y
   (for/list ([i (range n)])
     (observe (sample (list-ref y i)) (list-ref ys i))
     )
   
   ;;; (define post0 (normal mu0(0) 1))
  
   ; (list cp mu1 mu2)
   cp
   )
  )


(show-model (model *ys*) #:cred-mass 0.94)