Edit page
;;;fold: zip
(define (zip xs1 xs2) 
  (if (or (is_null xs1) (is_null xs2)) '() 
    (pair 
      (pair (first xs1) (pair (first xs2) '()))
      (zip (rest xs1) (rest xs2)))))
;;;

;; Matrix utility functions

(define (split-list-first split lst n)
  (if (= n 0)
      (cons split lst)
      (split-list-first (append split (list (car lst))) (cdr lst) (- n 1))))

(define (zero-vector n)
  ; list of n zeros
  (make-list n 0) )

;; Basic matrix functions

(define (m-reshape-col elements m n)
  ; Appends elements to matrix
  ; m rows, n columns
  ; Produces matrix as list of columns
  ; Assumes elements are column ordered i.e. (1 1) (2 1) ...
  (let* ([split (split-list-first '() elements m)]
         [new-column (car split)]
         [remaining-elements (cdr split)] )
    (if (= n 1)
        (cons new-column '())
        (cons new-column (m-reshape-col remaining-elements m (- n 1))) ) ) )

(define (m-reshape-row elements m n)
  ; Appends elements to matrix
  ; m rows, n columns
  ; Produces matrix as list of rows
  ; Assumes elements are row ordered i.e. (1 1) (1 2) ...
  (let* ([split (split-list-first '() elements n)]
         [new-row (car split)]
         [remaining-elements (cdr split)] )
    (if (= m 1)
        (cons new-row '())
        (cons new-row (m-reshape-row remaining-elements (- m 1) n)) ) ) )

(define (m->el A i j)
  ; Returns element (i j) of matrix A
  (list-ref (list-ref A (- i 1)) (- j 1)) )

(define (mr*v mr v)
  ; Matrix in row format right multiplied by vector
  (let ([dot-prod (sum (map * (car mr) v))])
    (if (= (length mr) 1)
        (list dot-prod)
        (cons dot-prod (mr*v (cdr mr) v)) ) ) )

(define (matrix-map-1 proc A)
  ; Map proc onto each element of A
  (map (lambda (a) (map proc a)) A) )

(define (matrix-map-2 proc A B)
  ; Map proc onto each element of A and B
  (map (lambda (a b) (map proc a b)) A B) )

;; Cholesky decomposition

;; A crude and slow implementation of the Cholesky-Banachiewicz algorithm

(define (chol A)
  ; Lower triangular Cholesky decomposition of A
  ; Assumes that A is a matrix of rows - but it might not matter
  (let ([L11 (sqrt (car (car A)))])
    (chol-square (chol-iter A (cons 2 1) (list (list L11)))) ) )

(define (chol-iter A ip L)
  ; Top level iteration of cholesky
  (if (> (car ip) (length A))
      L
      (chol-iter A (chol-next-index-pair ip) (chol-add-element A L (car ip) (cdr ip))) ) )

(define (chol-add-element P L i j)
  ; if j = 1 then add the new element into a new row
  ; else add the new element at the end of the last row
  (let ([new-el (chol-get-element P L i j)]
        [previous-rows (chol-truncate L)]
        [last-row (m->last-row/col L)])
    (cond
     ((= j 1) (append L (list (list new-el))))
     (else (append previous-rows (list (append last-row (list new-el))))) ) ) )

(define (chol-next-index-pair ip)
  ; Defines the order of iteration in Cholesky algorithm
  ; ip = (i . j)
  ; if i = j then return (i+1 . 1)
  ; else return (i . j+1)
  (cond
   ((= (car ip) (cdr ip)) (cons (+ (car ip) 1) 1))
   (else (cons (car ip) (+ (cdr ip) 1))) ) )

(define (chol-get-element P L i j)
  ; The Cholesky formula
  (cond
   ((= i j) (sqrt (- (m->el P i i) (chol-ps L i i))))
   (else (/ (- (m->el P i j) (chol-ps L i j)) (m->el L j j))) ) )

(define (chol-ps L i j)
  ; The sum in the formula for L_{ij}
  ; \sum_{k=1}^{j-1} L_{ik}L_{jk}
  (chol-ps-iter (- j 1) L i j) )

(define (chol-ps-iter k L i j)
  ; The main iterative loop for the Cholesky partial sum
  (cond
   ((= k 0) 0)
   ((= k 1) (* (m->el L i 1) (m->el L j 1)))
   (else (+ (* (m->el L i k) (m->el L j k))
            (chol-ps-iter (- k 1) L i j) )) ) )

;; Cholesky utilities

(define (chol-square A)
  ; pad the rows of A with zeros to make square
  (map (lambda (v) (append v (zero-vector (- (length A) (length v))))) A) )

;; GPs

(define (sqr x)
  (expt x 2) )

(define (mv-gaussian mean sigma)
  ; mean + chol(sigma) * gaussians
  (map + mean (mr*v (chol sigma) (repeat (length mean) (lambda () (gaussian 0 1))))) )

(define (euclidean-squared-distance X Y)
  (sum (map (lambda (x y) (sqr (- x y))) X Y)) )

(define (euclidean-squared-distances X Ys)
  (map (lambda (Y) (euclidean-squared-distance X Y)) Ys) )

(define (cov-se-iso-row x xs ell sf)
  (map (lambda (a) (* (exp (- 0 a)) (sqr sf))) 
       (map (lambda (b) (/ b (* 2 (sqr ell)))) 
            (euclidean-squared-distances x xs))) )

(define (cov-se-iso x ell sf)
  (map (lambda (a) (cov-se-iso-row a x ell sf)) x) )

(define (gp-se-iso x mean ell sf)
  (mv-gaussian mean (cov-se-iso x ell sf)) )

(define (sigmoid x)
  (/ (exp x) (+ 1 (exp x))) )

(define (m->last-row/col A)
  ; Returns the last row or column of a matrix
  ; Last row or column depends on the format of the matrix
  (car (reverse A)) )

(define (chol-truncate A)
  ; Returns everything but the last row
  (if (null? (cdr A)) 
      '()
      (reverse (cdr (reverse A))) ) )

(define input 
  '((1) (2) (3) (4) (5) (6) (7) (8) (9) (10)  (11) 
    (12) (13) (14) (15) (16) (17) (18) (19) (20)))

(define observations 
  '(#t #t #t #t #t #f #f #f #f #f #t #t #t #t #t #f #f #f #f #f))

(define samples
  (mh-query 
   
   10 20
   
   ;; Define the generative model and the likelihood of observation
   
   ;; Setup
   
   (define N 20)
   (define mu (zero-vector N))
   
   ; Broad priors on hyperparameters
   (define ell (+ 1 (min 2 (exp 0.1))))
   (define sf (+ 1 (exp 0.1)))
   
   ;; The model
   
   (define latent-gp (gp-se-iso input mu ell sf))
   
   ;; Likelihood calculations
   
   (define (bernoulli-likelihood p b)
     (if b p (- 1 p)) )
   
   (define (gp-class-log-lik gp obs)
     (sum (map (lambda (l b) (log (bernoulli-likelihood (sigmoid l) b))) gp obs)) )
   
   (define log-likelihood (gp-class-log-lik latent-gp observations))
   
   (define _ (factor log-likelihood))
   
   ;; Quantity to sample
   
   latent-gp
   
   #t
   
   ); mh-query
  ); define samples

(lineplot
 (zip (map car input) 
      (car (reverse samples))))

References: