# Probabilistic programming in POSIX C and OCaml: shift vs fork

Using probabilistic inference as an example, we liken the delimited control operator `shift` and the operating system call `fork`. Any POSIX programming language can thus be used for probabilistic modeling. In particular, we demonstrate probabilistic models and inference in ANSI C. This is joint work with Chung-chieh Shan.

## Probabilistic modeling in a general-purpose programming language

Our subject is developing models of the world using probability to account for inherent uncertainty, imperfect knowledge, inaccuracies in input data, or the abstraction over insignificant details. There is a lot of software for writing models, in graphical or textual forms, and performing sampling, parameter learning, inferring (unconditional) probability distribution of the outcome, computing outcome probabilities conditional on observations. In this article we concentrate on writing models as programs in a general-purpose programming language such as C, OCaml or Haskell. Probabilistic inference is then tantamount to running the program.

The canonical example is modeling a front lawn, whose grass may be wet because it rained or because the sprinkler was on. Since our knowledge is not perfect, we allow for 10% chance that the grass may be wet for some other reason. Also, if the rain was light, the grass might have dried up by the time we observe it; a sprinkler may too leave the grass dry if, say, water pressure was low. We admit these possibilities, assigning them probabilities 10% and 20%, respectively. Suppose we observe the grass is wet. What are the chances it rained?

The following is the model written as an OCaml program:

```     let grass_model () =   (* unit -> bool *)
let rain = flip 0.3 and sprinkler = flip 0.5 in
let grass_is_wet =
(flip 0.9 && rain) || (flip 0.8 && sprinkler) || flip 0.1 in
if not grass_is_wet then fail ();
rain
```
In the model, `rain` and `sprinkler` are independent boolean random variables, with discrete prior distributions `[(0.3, true); (0.7, false)]` for `rain` and `[(0.5, true); (0.5, false)]` for `sprinkler`. The variable `grass_is_wet` is the dependent random variable. The conditional statement asserts the evidence, using `fail ()` to report the contradiction of model's prediction with the observation. In the OCaml representation, `grass_model` is an OCaml function of the inferred type `unit->bool`; `rain`, `sprinkler` and `grass_is_wet` are regular OCaml variables; `flip` and `fail` are regular functions implemented below. We run the model as follows, obtaining the unnormalized posterior probability distribution of raining.
```     run_model grass_model;;
- : bool distrib = [(0.322, false); (0.2838, true)]
```
The sum of probabilities, 0.6058, is the `probability of evidence,' of the grass being wet. The answer to our question, the posterior probability of raining, is approximately `7/15`.
References
HANSEI: OCaml-embedded domain-specific language for probabilistic models and (nested) inference

Joyce, James: Bayes' Theorem
The Stanford Encyclopedia of Philosophy, Edward N. Zalta (ed.)
< http://plato.stanford.edu/entries/bayes-theorem >

## The lawn model in C

We can write the same lawn model as a C program:
```     bool grass_model (void) {
bool rain      = flip(0.3);
bool sprinkler = flip(0.5);
bool grass_is_wet =
(flip(0.9) && rain) || (flip(0.8) && sprinkler) || flip(0.1);
if( !grass_is_wet )
fail();
return rain;
}
```
where the functions `flip` and `fail` are defined below. Modulo minor syntactic differences, the code looks quite like the one in OCaml. We run the model using this function:
```     void runit(void) {
bool result = grass_model();
printf("%g %c\n",MyWeight,result ? 'T' : 'F');
}
```
The standard output of the program is a sequence of lines consisting of two fields separated by a space. The first field is the probability, the second field is the outcome, the letter T or F. The output needs to be post-processed and collated -- for example, by piping the output to the following AWK program. In the spirit of UNIX, each program ought to do one thing: generate results or post-process them.
```     awk '{r[\$2] += \$1}; END {printf "T: %g; F: %g",r["T"],r["F"]}'
```
The produced output is: `T: 0.2838; F: 0.322`
Version
The current version is 1.1, Jan 29, 2009.
References
grass.c [3K]
The complete C source code

## Probabilistic programming in OCaml

To run the OCaml program representing the lawn model, we need to implement the functions `flip`, `fail` and `run_model`. The function `flip` implements a flip of a weighted coin. We may also regard `(flip p)` to be a random boolean value, assuming the value true with the probability `p`. Finally, we may treat a program as representing a distribution, a weighted list of possible outcomes. A deterministic program with the outcome `v` denotes the distribution `[(1.,v)]`. Then `(flip p)` represents the distribution `[(p,true); (1-p,false)]`. The function `fail` is too non-deterministic: it denotes the impossibility, the distribution `[]`. The function `run_model` takes a program as the argument and converts, or reifies, it to the probability distribution the program represents.

Here is the complete code, which relies on the library delimcc of delimited continuations in OCaml.

```     open Delimcc

type prob = float
type 'a distrib = (prob * 'a) list

let scale_prob pscale : 'a distrib -> 'a distrib =
List.map (fun (p,x) -> (p *. pscale,x))

let collate : 'a distrib -> 'a distrib =
let collator collated (p,v) =
match List.partition (fun (_,v') -> v = v') collated with
| ([],      l) -> (p,      v)::l
| ([(p',_)],l) -> (p +. p',v)::l in
List.fold_left collator []

let p0 = new_prompt ()

let flip p =  (* float -> bool *)
shift p0 (fun k -> scale_prob p (k true) @
scale_prob (1. -. p) (k false))
let fail () = abort p0 []

let run_model model =
collate (push_prompt p0 (fun () -> [(1., model ())]))
```
The type annotations are all optional but given for clarity. It is instructive to compare the above code with the implementation in C below: whereas `flip` in OCaml uses `shift`, `flip` in C uses `fork(2)`. In both languages, `fail` is just abort.
References
Delimited continuations in OCaml

## C as a probabilistic programming language

The idea of a program representing a distribution of its possible outcomes can also be implemented in C. The outcome of a deterministic program can be computed by one process. For a program with several possible outcomes we can use several processes to compute the outcomes in parallel. Therefore, to evaluate `(flip p)` we should fork the current process in two. We track the weight of each outcome in a per-process global variable `MyWeight`. The function `fail` is just the synonym for `abort(3)`.
```     static double MyWeight = 1.0;   /* changed at the start-up of each process */

bool flip(const double p) {
pid_t pid = fork();
if(pid == 0) {                /* in the child process, to handle False */
MyWeight *= 1.0 - p;
return false;
}
else if(pid == (pid_t)(-1)) {
perror("fork failed");
abort();
} else {                      /* in the parent process, to handle True */
MyWeight *= p;
return true;
}
}

int main(void) {
int status;
runit();
while (wait(&status) > 0)
;
return 0;
}
```
Version
The current version is 1.1, Jan 29, 2009.
References
grass.c [3K]
The complete C source code

## Discussion

Every Unix programmer already knows what a delimited continuation is: it is the `state' of an unfinished user process, captured on interrupt or system call, stored in a process table, written on disk upon a checkpoint, sent across the network during the process migration, discarded on abort, or resumed by the scheduler, sometimes twice (when handling `fork`). Conversely, if a programming language provides for capturing delimited continuations, we can implement our own OS with our own multi-processing. The ZFS project did exactly that.

Our library for probabilistic programming in OCaml is richer. It lets us evaluate the same model in different ways, obtaining either the exact distribution of model's outcomes, a sample outcome, or a sequence of sample outcomes, from which we approximate their distribution. The users of our library may write their own optimized inference and sampling strategies, necessary to deal with realistic models. In addition to the basic `flip` and `fail`, our library provides the third primitive, to reify a probabilistic program into a lazy search tree, which we can then incrementally explore. On a POSIX system, we may implement such a primitive by using, for example, a programmable debugger or by writing a system call interceptor.

References
Delimited continuations in operating systems

Zipper-based file server/OS

HANSEI: OCaml-embedded domain-specific language for probabilistic models and (nested) inference

### Last updated February 1, 2009

This site's top page is http://okmij.org/ftp/

oleg-at-okmij.org