An unknown-order polynomial is searched for here, one that plausibly generated a small set of x-y data pairs. A random order between one and four gets drawn along with matching coefficients, and Metropolis-Hastings conditions on the data to infer the polynomial most likely to have produced it.
(define (zip xs1 xs2)
(if (or (null? xs1) (null? xs2)) '()
(pair
(pair (first xs1) (pair (first xs2) '()))
(zip (cdr xs1) (cdr xs2)))))
(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 samples
(mh-query 10 100
(define (randint l h) (+ l (sample-integer (+ (- h l) 1))))
(define order (randint 1 4))
(define (my-range n) (if (equal? n 0) '() (pair (- n 1) (my-range (- n 1)))))
(define (my-range2 n) (reverse (my-range n)))
(define (pow n x) (if (equal? n 0) 1.0 (* x (pow (- n 1) x))))
(define expts (my-range2 order))
(define coefs (repeat order (lambda () (uniform -3.0 3.0))))
(define data
'((-4.0 16.0) (-3.0 9.0) (-2.0 4.0) (-1.0 1.0)
(0.0 0.0) (1.0 1.0) (2.0 4.0) (3.0 9.0) (4.0 16.0)))
(define (poly-fx x)
(sum (map (lambda (c-n) (* (first c-n) (pow (second c-n) x)))
(zip coefs expts))))
(define soft-eq (lambda (x y) (factor (glp x 0.1 y))))
(define pairfactors (map (lambda (x-y) (soft-eq (second x-y) (poly-fx (first x-y)))) data))
(define sample (pair order coefs))
sample
#t))
samples
References: