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.
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.Joyce, James: Bayes' Theorem
The Stanford Encyclopedia of Philosophy, Edward N. Zalta (ed.)
< http://plato.stanford.edu/entries/bayes-theorem >
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.322flip, 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.(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;
}
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.
HANSEI: OCaml-embedded domain-specific language for probabilistic models and (nested) inference
oleg-at-okmij.org
Your comments, problem reports, questions are very welcome!
Converted from HSXML by HSXML->HTML