Edit page
Note: The code on this page may fail. (The MH chain does not mix.)

Bayesian inference can be used to simultaneously find the order of a polynomial and its coefficients, trading off fit against prior preference for simpler models.

;;;fold: my-pi, glp, range
(define my-pi 3.14159265358979323)

(define (glp mean ssq smp)
  (let ([diff (- smp mean)])
    (- (- (/ (* diff diff) (* 2.0 ssq)))
       (* 0.5 (+ (log 2) (log my-pi) (log ssq))))))

(define (xrange i j)
  (if (= i j)
      (cons i (range (+ i 1) j))))

(define (range i j)
  (reverse (xrange i j)))
(define (make-poly c)
  (lambda (x) (apply + (map (lambda (a b) (* a (expt x b)))
                            (iota (length c))))))

(define x-vals (map (lambda (x) (/ x 1)) (range -5 5)))

(define true-coeffs '(-.5 1 .5))

(define true-y-vals (map (make-poly true-coeffs) x-vals))

(define obs-noise 0.2)

(define obs-y-vals
  (map (lambda (x) (gaussian x obs-noise)) true-y-vals))

(define (noisy-equals? x y)
  (factor (glp x obs-noise y))

(define samples
   10000 5
   (define poly-order (sample-integer 4))
   (define coefficients
     (repeat (+ poly-order 1)
             (lambda () (gaussian 0 2))))
   (define y-vals
     (map (make-poly coefficients) x-vals))
   (list poly-order coefficients)
   (all (map noisy-equals? y-vals obs-y-vals))))

(hist (map first samples))