ece59500nl/probabilistic-scheme.ss at main · qobi/ece59500nl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
(define first car)

(define rest cdr)

(define (flip alpha)
 (unless (<= 0 alpha 1) (error "Alpha not probability"))
 (cond ((= alpha 0) #f)
       ((= alpha 1) #t)
       (else (call-with-current-continuation
	      (lambda (c)
	       (let ((saved-bottom bottom))
		(set! bottom
		      (lambda ()
		       (set! bottom saved-bottom)
		       (c #f)))
		#t))))))

(define (bottom) (error "Top-level bottom"))

(define (sample x)
 (let loop ((x (normalize x)) (p 1))
  (cond ((or (null? x) (zero? p)) (bottom))
	((flip (/ (cdr (first x)) p)) (car (first x)))
	(else (loop (rest x) (- p (cdr (first x))))))))

(define (remove-if-not p l)
 (let loop ((l l) (c '()))
  (cond ((null? l) (reverse c))
	((p (first l)) (loop (rest l) (cons (first l) c)))
	(else (loop (rest l) c)))))

(define (map-reduce g i f l)
 (if (null? l) i (map-reduce g (g i (f (car l))) f (cdr l))))

(define (normalize distribution)
 ;; It doesn't matter that we don't coalesce since except for efficiency
 ;; you can say a distribution with duplicates in its support is just an
 ;; implicit representation of a proper distribution.
 (let ((distribution
	(remove-if-not (lambda (pair) (> (cdr pair) 0)) distribution))
       (sum (map-reduce + 0 cdr distribution)))
  (map (lambda (pair) (cons (car pair) (/ (cdr pair) sum))) distribution)))

(define-syntax distribution
 (syntax-rules ()
  ((distributions e)
   (normalize
    (call-with-current-continuation
     (lambda (c)
      (let ((values '())
            (saved-flip flip)
            (saved-bottom bottom)
            (p 1))
       (set! flip
             (lambda (alpha)
              (unless (<= 0 alpha 1) (error "Alpha not probability"))
	      (cond ((= alpha 0) #f)
		    ((= alpha 1) #t)
		    (else (call-with-current-continuation
			   (lambda (c)
			    (let ((saved-bottom bottom)
				  (saved-p p))
			     (set! bottom
				   (lambda ()
				    (set! bottom saved-bottom)
				    (set! p (* (- 1 alpha) saved-p))
				    (c #f)))
			     (set! p (* alpha p))
			     #t)))))))
       (set! bottom
             (lambda ()
              (set! flip saved-flip)
              (set! bottom saved-bottom)
              (c (reverse values))))
       (set! values (cons (cons e p) values))
       (bottom))))))))

(define-syntax probability
 ;; marginal probability
 (syntax-rules ()
  ((probability e) (map-reduce + 0 cdr (remove-if-not car (distribution e))))))

(define-syntax support
 (syntax-rules () ((support e) (map car (distribution e)))))

(define-syntax expected-value
 ;; expectation
 (syntax-rules ()
  ((expected-value e)
   (map-reduce +
	       0
	       (lambda (pair) (* (car pair) (cdr pair)))
	       (distribution e)))))

(define-syntax variance
 (syntax-rules ()
  ((variance e) (- (expected-value (expt e 2)) (expt (expected-value e) 2)))))

(define-syntax raw-moment
 (syntax-rules () ((raw-moment n e) (expected-value (expt e n)))))

(define-syntax central-moment
 (syntax-rules ()
  ((central-moment n e)
   (let ((mu (expected-value e))) (expected-value (expt (- e mu) n))))))

(define (lg x) (/ (log x) (log 2)))

(define-syntax entropy
 (syntax-rules ()
  ((entropy e)
   (- (map-reduce +
		  0
		  (lambda (pair) (* (cdr pair) (lg (cdr pair))))
		  (distribution e))))))

(define-syntax kl-divergence
 ;; Kullback-Leibler divergence
 ;;\needswork: abstract equal?
 ;;            check that (set-equal? equal? (support e1) (support e2))
 (syntax-rules ()
  ((kl-divergence e1 e2)
   (map-reduce +
	       0
	       (lambda (x)
		(let ((p (probability (equal? e1 x)))
		      (q (probability (equal? e2 x))))
		 (* p (lg (/ p q)))))
	       (support e1)))))

(define (product-distribution f d1 d2)
 (distribution (f (sample d1) (sample d2))))

(define (marginal-distribution f d)
 (map (lambda (pair) (cons (f (car pair)) (cdr pair))) d))

(define-syntax mutual-information
 ;;\needswork: abstract car cdr cons
 (syntax-rules ()
  ((mutual-information e)
   (let* ((joint-distribution (distribution e))
	  (marginal-distribution-car
	   (marginal-distribution car joint-distribution))
	  (marginal-distribution-cdr
	   (marginal-distribution car joint-distribution)))
    (kl-divergence
     e
     (sample (product-distribution cons
				   marginal-distribution-car
				   marginal-distribution-cdr)))))))