sugar observe from, dockerfile, readme

README.md

@@ -1,12 +1,12 @@
 # Simple Imperative Minimal Probabilistic Language
-This is a toy probabilistic programming system written in beautiful racket.
-The goal is to provide an off-the-shelf imperative programming language with built-in language constructs for statistical computations such as sampling and bayesian inversion.
+This is a toy probabilistic programming system written in beautiful racket, which is by no means efficient or even safe.
+The goal is to gain more insight in the mechanics of step-by-step probabilistic computation, with language constructs for sampling and bayesian inversion.
 
 ## Installation Procedures
 Install docker on your machine if not running Linux.
-The system has only been tested on Debian (bullseye) and Arch Linux.
+The system has been tested on Debian (bullseye) and Arch Linux.
 
-### In Linux
+### On Linux
 Make sure you have installed `git`, `racket`, and `sqlite3`.
 Also install `beautiful-racket` by executing the following command:
 ```
@@ -16,39 +16,34 @@ After cloning this repository and changing the working directory to it, run
 ```
 $ raco pkg install
 ```
-to install simpl.
-Run the unit (some are probabilistic and might give errors) tests with
-```
-$ raco test tests
-```
+to install the programming system.
 
 ### Using Docker
-After cloning this repository, build a docker image by executing the following command within the repository:
+After cloning this repository, build a docker image by executing the following command within the local repository:
 ```
-$ docker build . 
+$ docker build -t debian:simpl . 
 ```
-
-If successful, inspect the IDs of all the images you have to obtain your appropriate image ID `N`.
-This can be done by executing
+If successful, the docker image you just installed shows up in the list of images:
 ```
 $ docker image list
 REPOSITORY    TAG     IMAGE ID    ...
-<none>        <none>  N
+debian        simpl   N           ...
 ```
-Use this ID `N` to execute the image as follows (only the first couple of digits is required to enter):
+You can run the image in an interactive container by doing
 ```
 $ docker run -it N
 ```
 
 ## Testing
-Run the unit tests with the command below.
-NB: some tests fail with a probability of approximately 5%; most test runs should succeed though, so rerunning the tests after a failure likely "fixes" things.
+Run the unit tests with 
 ```
 $ raco test tests
 ```
+Note that in some test cases, sampled values are checked to be within a certain range, and may fail with some low but nonzero probability.
+This is to be expected.
 
 ## Examples
-There is some examples available in the directory `examples`. 
+There are some examples available in the directory `examples`. 
 
 ### Normal sampling
 The file `normal-sampling.rkt`, for example, samples 1000 times from a normal(100,20) distribution and calculates the sample mean and sample variance.  
@@ -57,10 +52,23 @@ $ racket examples/normal-sampling.rkt
 ```
 The result should be _approximately_ (100 20).
 
+The first line of `normal-sampling.rkt` is 
+```
+#lang simpl
+```
+This tells racket that it is supposed to interpret the rest of the file as the simple imperative minimalist programming language that is this tool.
+
 ### Burglary alarm
-The file `burglary.rkt` contains a probabilistic model for a simple Bayesian inference procedure on Bernoulli distributed random variables representing the occurence of a burglary and an alarm going off.
+The file `gender-height.rkt` contains a probabilistic model for a simple Bayesian inference procedure on the Bernoulli distributed random variable for the gender of a person, conditioned on the fact that we observe a height of 190 cm.
+```
+$ racket examples/gender-height.rkt
+```
+The first line contains `#lang simpl/sybex`, which tells racket to do symbolic execution, rather than concrete. 
+The result will thus (always) be
 ```
-$ racket examples/burglary.rkt
+(list
+ (end-config 1 9.684491216181099e-5 '((< X_0 0.496)))
+ (end-config 0 0.008598284478336488 '((not (< X_0 0.496)))))
 ```
-Out of 10000 executions, the result will 9999 times be (0.001 0), meaning that outcome 0 under this model has a weight (unnormalized probability) of 0.001.
-Very rarely, the result is (0.95 1), meaning that outcome 1 has a weight (unnormalized probability) of 0.95.
+The structures `end-config` contain respectively the _return value_ (1 and 0 here), the _likelihood_ of the outcome, and the path condition.
+The path condition is expressed using symbolic random variables, in this case X_0.

dockerfile

@@ -1,4 +1,7 @@
 FROM debian:bullseye
+LABEL maintainer="erikvoogd@ifi.uio.no"
+LABEL build_date="2022-06-30"
+
 WORKDIR /code
 RUN apt-get update 
 RUN apt-get install -y git racket sqlite3

examples/burglary.rkt

@@ -5,5 +5,5 @@ if (b == 1) {
 } else {
   a := 0.001
 };
-score (a);
+observe 1 from bern(a); // sugar for `score(a)`
 return b

examples/gender-height.rkt

@@ -0,0 +1,11 @@
+#lang simpl/symbex
+gender ~ bern(0.496);
+if (gender == 0) {
+  mean := 175;
+  var  := 64
+} else {
+  mean := 165;
+  var  := 49
+};
+observe 190 from normal (mean,var);
+return gender

examples/phonecalls.rkt

@@ -13,7 +13,6 @@ if (x == 1) {
   r:=3          // we expect 3 calls
 };
 
-// we observe 4 from poisson(r) (use density function)
-score (r^4*2.718281828459045^(-r)/24); 
+observe 4 from poisson(r);
 
 return x

expander.rkt

@@ -36,6 +36,8 @@
        `(while ,(expand-expr c) ,(expand-stmts body))]
       [(list 'stmt "score" "(" e ")")
        `(score ,(expand-expr e))]
+      [(list 'stmt "observe" e1 "from" d e2)
+       `(score ,(score-sugar d (expand-expr e1) (expand-expr e2)))]
       [(list 'stmt "return" e)
        (list 'return (expand-expr e))])))
 
@@ -144,3 +146,27 @@
       [(? symbol? x) x]
       [(? boolean? x) x]
       [(? string? x) x])))
+
+(define score-sugar
+  (lambda (d e1 e2)
+    (match d
+      ("bern"
+       (if (list? e2)
+           (error (format "observe statement with bernoulli should have one parameter: `~a`" e2))
+           `(* (expt ,e2 ,e1) (expt (- 1 ,e2) (- 1 ,e1)))))
+      ("poisson"
+       (if (list? e2)
+           (error (format "observe statement with exponential distribution requires one parameter: `~a`" e2))
+           (let ((e 2.71828182845904523536))
+             `(* (/ (expt ,e2 ,e1) (fac ,e1)) (expt ,e (- ,e2))))))
+      ("normal"
+       (if (not (and (list? e2) (eq? (length e2) 2)))
+           (error (format "observe statement with normal should have two parameters: `~a`" e2))
+           (let ((mu (car e2))
+                 (var (cadr e2))
+                 (sqr (lambda (x) (* x x)))
+                 (e 2.71828182845904523536)
+                 (pi 3.14159265358979323846))
+             `(/ (expt ,e (- (/ (* (- ,e1 ,mu) (- ,e1 ,mu)) (* 2 ,var)))) (sqrt (* (* 2 ,pi) ,var))))))
+      (else
+       (error (format "observe from ~a not implemented" d))))))

lexer.rkt

@@ -12,7 +12,7 @@
     [(:or (from/stop-before "//" "\n")
           (from/to "/*" "*/"))
      (token lexeme #:skip? #t)]
-    [(:or "skip" "abort" "if" "else" "while" "true" "false" "return" "array" "score")
+    [(:or "skip" "abort" "if" "else" "while" "true" "false" "return" "array" "score" "observe" "from")
      (token lexeme lexeme)]
     [(:or "uniform" "bern" "binom" "poisson" "normal" "exp")
      (token 'DIST lexeme)]

parser.rkt

@@ -10,6 +10,7 @@ stmt : "skip"
      | "while" "(" be ")" "{" stmts "}"
      | "score" "(" e ")"
      | "return" e
+     | "observe" e "from" DIST e
 
 be  : be1 [ "||" be ]
 be1 : be2 [ "&&" be1 ]

sym-substeval.rkt

@@ -20,7 +20,9 @@
       [(list (? bin-op? op) left right)
        (list op (subst sub left) (subst sub right))]
       [(list (? un-op? op) opnd)
-       (eval (list op (subst sub opnd)) ns)]
+       (if (eq? op 'fac)
+           (fac (subst sub opnd))
+           (list op (subst sub opnd)))]
       [(? list? x) (map (lambda (y) (subst sub y)) x)]
       [(? number? x) x]
       [(? boolean? x) x]

utils.rkt

@@ -3,7 +3,11 @@
 (define add-op? (lambda (x) (memq x '(+ -))))
 (define mul-op? (lambda (x) (memq x '(* /))))
 (define bin-op? (lambda (x) (memq x '(+ - * / || && < <= > >= equal? expt))) )
-(define un-op?  (lambda (x) (memq x '(not - sqrt))))
+(define un-op?  (lambda (x) (memq x '(not - sqrt fac))))
+
+(define fac
+  (lambda (x)
+    (apply * (build-list x add1))))
 
 (define (round-off z n)
   (let ((power (expt 10 n)))
@@ -30,5 +34,5 @@
       [(? symbol? x) #t])))
 
 (provide add-op? mul-op? bin-op? un-op? 
-         round-off 
+         round-off fac
          make-symbolic symbolic-variable? symbolic-expr?)