<jupyter_start><jupyter_text>Getting started with probabilistic programming in GambleThis document introduces the basic ideas of probabilistic programming using Gamble, andprovides lots of small examples as starting points for your own projects.  You can tweak and re-run all the code blocks.This primary author of this introduction is <NAME> (<EMAIL>), who welcomes bug reports and suggestions for improvement. PrerequisitesThis material assumes that you are familiar with elementary probability theory, that youknow how to program in some language, and that you are at least a little bitfamiliar with the Scheme programming language.  If your probability isrusty, there are nice reviews [here](http://cs229.stanford.edu/section/cs229-prob.pdf) and[here](http://www.math.wisc.edu/~anderson/605F11/Notes/StochBioChapter2.pdf).For a more leisurely discussion, with a focus on the Bayesian perspective, and some modeling examples, click[here](http://cocosci.berkeley.edu/tom/papers/bayeschapter.pdf). For help with Scheme, and, in particular,the Racket dialect spoken here, [this](http://docs.racket-lang.org/guide/index.html) is a good place to start.[TOC]  Contents* [What is probabilistic programming and why do we need it?](What is probabilistic programming and why do we need it?)* [Inference algorithms](Inference algorithms)* [Beyond coin flips](Beyond coin flips)* [Indexed sets of random variables](Indexed sets of random variables)  What is probabilistic programming and why do we need it?What is probabilistic programming?Computer programs are a familiar and flexible way to represent processes, and running those programs repeatedly is a straightforward way to simulate those processes. Like any other formal mathematical *model*, a model in program form abstracts away (currently) unimportant details so that we can describe the essentials of a system, and derive properties that were not obvious without the model. A computer program is a model with the additional virtue of being *executable*.  A *probabilistic* (or *stochastic*) model is a model that includes random elements. Probabilistic models are incredibly useful, even though truly probablistic systems are rare (even coin flips and dice rolls are deterministic, given enough information). The probabilistic parts are summaries of our ignorance of the many factors that we would have to know precisely to model (e.g.) coin flipping as deterministic.Given thatall popular programming languages include procedures for drawing (pseudo-)randomnumbers, programs can easily represent random processes, including complex, multi-stepprocesses with both random and deterministic parts.  We will call programs with random expressions "probabilistic programs" (PPs, for short).  So,>**Probabilistic programs are programs with random components.** (If this was all there was to it, of course, we wouldn't need new languages.  We'll get to why we need something new in the next section.)We can think of a complete program as a procedure (often one that takes no arguments).  Running this programmany times will produce a sample resembling the probability distribution over executionsdefined by the program. Before we start...When you first open a new notebook, it's a *Racket* notebook.  Racket is the language (or language-building toolkit) in which Gamble is defined.  To make it possible to run Gamble programs, the first thing you will want to do is this:<jupyter_code>(require gamble)<jupyter_output><empty_output><jupyter_text>If you don't evaluate a (require gamble) expression first, you won't be able to evaluate any Gamble code. The first Gamble programWith that bit of plumbing taken care of, consider the simplestprobabilistic "process": a single coin toss.  We can model it with a random variable thattakes on the values $\verb+t+$ and $\verb+f+$ (Racket's Boolean values) with equal probability.  We represent this "system" with the smallest possible Gamble program:<jupyter_code>(flip)<jupyter_output><empty_output><jupyter_text>If we run this program many times, it will return $\verb+t+$ about half the time and $\verb+f+$ the other half.  The program below flips the coin100 times and counts the $\verb+t+$s:<jupyter_code>(length 
 (filter (lambda (x) x)
         (repeat flip 100)))<jupyter_output><empty_output><jupyter_text>That's pretty much the "Hello, world!" of probabilistic programming.  Here's the line-by-line breakdown: repeat executes the **thunk** (that's a procedure that takes no arguments) passed in as its first argument $n$ times, where $n$ is an integer value passed in as the second argument.  filter runs a predicate (first argument) over a list (second argument), and returns just those elements of the list for which the predicate evaluates to t.  Even if you haven't seen it before, you can probably guess what length does.flip is one of about 20 primitive distibutions, which can be composed with each other and with other (pure) Racket constructs to make Gamble models.  We'll get to a few of the most common ones in this introduction, and the rest are documented here (ref).There is another primitive distribution very much like flip, except that it yields 0 or 1 instead of t or f.<jupyter_code>(bernoulli)<jupyter_output><empty_output><jupyter_text>flip is useful for random variables used as conditions in if expressions, but counts and averages of many "flips" are easier with bernoulli.<jupyter_code>(apply + (repeat bernoulli 100))<jupyter_output><empty_output><jupyter_text>Conditional distributionsBigger computer programs are composed of smaller ones.  So are Gamble programs representing complex joint probability distributions.  flip and bernoulli, above, are passed no arguments, and default to equal probability for both outcomes.  If all variables in a joint distribution are produced by procedures of no arguments, there are no attachement points between variables, so no interesting joint distribution can be constructed.  If we think in terms of cause and effect (a risky but useful heuristic), the zero-argument version of flip is an isolated event, without cause.  It feels no pain.  It never cries.  The bread and butter of probabilistic programming is *conditional* distributions, which map *parameters* to distributions.For example, flip can take a bias argument to make a "trick coin".<jupyter_code>(flip 0.9)<jupyter_output><empty_output><jupyter_text>This is important, so I will say it again: when procedures with parameters are (or contain) random choices, like flip, they represent *conditional* distributions--that is, functions from arguments (here, a coin's bias towards t) todistributions (or samplers of distributions).If we flip such a coin many times, it behaves just as you'd expect.<jupyter_code>(length 
  (filter (lambda (x) x)
          (repeat (lambda () (flip 0.9)) 100)))<jupyter_output><empty_output><jupyter_text>**Syntax aside:** the expression is uglier than it was for zero-argument version of flip, because repeat requires as its first argument a thunk that can be evaluated to produce a sample.  (flip 0.9), withoutthe enclosing lambda, is not a thunk, but, rather, just a single value--either t or f.  We could use a few defines to make this less ugly:<jupyter_code>(define (flip-mostly-heads-coin) (flip 0.9))

(define (is-true? x) x)

(length (filter is-true? (repeat flip-mostly-heads-coin 100)))<jupyter_output><empty_output><jupyter_text>Gluing several random choices togetherAbove, we reached in and set the bias parameter ourselves, but the bias itself could be a random choice.   Here is a Gamble program that actually **links** several random variables, by sharing a bias variable coin-weight across the two flip expressions:<jupyter_code>(define (two-weighted-flips) 
  (let ([coin-weight (uniform 0 1)])
  (list (flip coin-weight) (flip coin-weight))))

(two-weighted-flips)<jupyter_output><empty_output><jupyter_text>two-weighted-flips doesn't look like much, but it's actually a big step.  It is through these linkages that information flows from one part of a mode to another.  Finding out the value of a coin flip, here, tells us something about the coin weight, and vice versa.  This process of incorporating observations of some variables to (rationally) change beliefs about other, *un*observed variables--*inference*--is where things gets tricky. Why we need something newSo far, we haven't done anything we couldn't have done in just about any existing programming language. But inference is where probabilistic programming in ordinary programming languages fails.  In short, **sampling from the possible executions of the process doesn'tsolve the kinds of problems people usually want to solve**. We usually **know something** about how the process executed, and want to know other things.*Some* inferences *are* easy.  For example we can easily sample many values of two-weighted-flips, and look at (say) the proportion of heads on the second flip, or the proportion of the time the two flips are equal.<jupyter_code>(require racket/list)

(length
 (filter (lambda (ht-list)
           (equal? (first ht-list)
               (second ht-list))) (repeat two-weighted-flips 100)))<jupyter_output><empty_output><jupyter_text>The (require racket/list) is so I can use first and second, which I prefer to the antique lisp forms car and cadr.We can also easily sample from the model in which a "root" random variable--in this case,the coin weight--has been observed.  This is just requires a bit of "program surgery", replacing the random draw of the weight with assignment to a fixed value:<jupyter_code>(define (two-mostly-heads-flips) 
  (let ([coin-weight 0.9])
  (list (flip coin-weight) (flip coin-weight))))

(two-mostly-heads-flips)<jupyter_output><empty_output><jupyter_text>Now let's consider a harder question: having observed the first flip's outcome as t, in what proportion of samples (of the whole process) does the other flip also have value t?  We can't just fix the value of the first flip, because that severs its connection to the coin weight.  In general, then, inference is hard, and we can't do it by simply evaluating a program written in a general-purpose programming language.First, we consider two simple inference strategies, both of which are implemented in Gamble.Inference strategy 1: reject samples that don't agree with the observationRejection sampling is dead simple, and often used as the basis for an operational semantics for probabilistic programming (ref).  The idea is to generate many samples of all the random variables, from the original model--the "prior", in Bayesian lingo--and simply throw out those for which the sampled value(s) of the observed random variable(s) don't agree with the observed value(s).  We could easily roll our own using filter, but rejection sampling is built in to Gamble, and provides a nice introduction to the general pattern of a query.Instead of editing the model, the normal Gamble inference workflow specifies a sampler object, defined by* a generative model (the program that specifies the process),* some observations, and* a query expression.The class of the object itself specifies the sampling method.<jupyter_code>(define two-wtd-cn-rej-sampler
  (rejection-sampler
   
    ;; Model.
    (define coin-weight (uniform 0 1))
    (define flip1 (flip coin-weight))
    (define flip2 (flip coin-weight))
   
    ;; Observations.
    (observe/fail flip1 #t)
   
    ;; Query ("return") expression.
    flip2))
   
   
(length (filter is-true? (repeat two-wtd-cn-rej-sampler 10000)))<jupyter_output><empty_output><jupyter_text>It looks like observing the first flip to be t tells us that the second is t about two thirds of the time.  As we will see in a bit, this happens to be one of the few special cases where we can verify that answer against a closed-form solution.TODO: Show more idiomatic, vector-based sample generating/handling procedures.Inference strategy 2: weight samples of other variables by the degree to which they "agree with" the observationRejection is a wasteful procedure, especially when the observed values occur in only a small fraction of samples from the prior.  "Likelihood weighting" is a bit smarter; instead of actually drawing values for the observed random variables, many or most of which will be incompatible with the observed values, it *weights* samples of the *other* random variables by how probable the observed values are, given those sampled values.<jupyter_code>(define two-wtd-cn-lw-sampler
  (importance-sampler
   
    ;; Model.
    (deflazy coin-weight (uniform 0 1))
    (deflazy flip1 (flip coin-weight))
    (deflazy flip2 (flip coin-weight))
   
    ;; Observations.
    (observe flip1 #t)
   
    ;; Query ("return") expression.
    flip2))


(generate-weighted-samples two-wtd-cn-lw-sampler 30)<jupyter_output><empty_output><jupyter_text>Indexed sets of random variablesIf we define a to have random value $v$, then subsequent references to a will have value $v$.  Very often, we will have a *set* of random variables, each indexed by some identifier--say, the heights of a number of people--that are drawn similarly.  Instead of giving each variable in this set a unique name and drawing a value for it with a define statement, we'd like to do something like a for loop over an array.  In scheme, a natural generalization of a simple random variable to an indexed set is a function. An array is just a function from integers to values, but we need not limit ourselves to integer indices.  We might try something like this:    (define (height person) (normal 69 3))    However, there is a problem: a function is evaluated each time it is encountered, so several references to (height "Bob") will, in general, have different random values.  This is not what we mean.  We want there to be just one random value for (height "Bob"), for all references.  What we want is the special Gamble syntax defmem, used like this:    (defmem (height person) (normal 69 3))    I like to think of this as drawing an infinite table assigning each possible input to a normal draw.  Subsequent references are just table lookups.  Of course, the implementation has to be clever and only draw values for table entries that are actually used, when they are first needed.  The "mem" in defmem is for "memoize", which just means wrapping a function so that only the first instance of a call (with a particular argument) is actually evaluated, with subsequent calls looking up the value generated the first time.  With deterministic procedures, memoization is a speed-up technique, but here it gets us from the wrong program meaning to the right one. TODO: make this work, and then answer a question using weighted samples. Also, successfully switch to deflazy and observe (not observe/fail), and explain the differences between the two "observes"...Maybe first do with observe/fail, then note that this is essentially rejection. Then introduce the notion of an observable context...observe versus observe/failobserve is better than observe/fail in just the same way that likelihood weighting is better than rejection sampling: rather than waiting to see what value has been drawn, to apply an observation, observe weights the rest of the draw by the observed value.  However, observe has some restrictions on its use.  First of all, it can't be used when a (possibly different) value for the observed random value has already been drawn!  By substituting deflazy for define, [somewhere], we tell the compiler[/interpreter?] not to draw that value immediately, but wait, instead, until it is required to produce the query value or score a state [...].  A deflazy and an observe of the same variable can be combined into a weighting expression by the compiler. [...Ryan...].Explore the space of values of unobserved random variables by a sequence of random *proposals*The workhorse of inference, for hard problems, is the *Metropolis-Hastings* sampler.  Why it works is beyond the scope of this introduction, but the intuition is not too hard.* The sampled and observed values define a *score* (probability or log-probability) for every state.* From any current state, simple random perturbations of values can *propose* a new state.* The score function can be used to decide whether to accept the proposed state or keep the old one.* As if by magic, this procedure eventually samples from the correct posterior distribution.The "eventually" part means that it is usually a good idea to sample for a while (a "burn-in" period) before collecting any values.  Also, because successive samples are far from independent, it can be space efficient to "thin" the sample sequence--that is, to skip some number of samples between recorded ones.  These two delays are input to generate-samples by the optional, keyword-identified arguments :burn and :thin.<jupyter_code>(require racket/vector)

(define two-wtd-cn-mh-sampler
  (mh-sampler
   
    ;; Model.
    (define coin-weight (uniform 0 1))
    (define flip1 (flip coin-weight))
    (define flip2 (flip coin-weight))
   
    ;; Observations.
    (observe/fail flip1 #t)
   
    ;; Query ("return") expression.
    flip2))


(time (vector-length 
 (vector-filter is-true? 
                (generate-samples two-wtd-cn-mh-sampler 100 #:burn 100 #:thin 100))))
(require racket/vector)

(define two-wtd-cn-mh-sampler2
  (mh-sampler
   
    ;; Model.
    (define coin-weight (uniform 0 1))
    (deflazy flip1 (flip coin-weight))
    (define flip2 (flip coin-weight))
   
    ;; Observations.
    (observe flip1 #t)
   
    ;; Query ("return") expression.
    flip2))


(time (vector-length 
 (vector-filter is-true? 
                (generate-samples two-wtd-cn-mh-sampler2 100 #:burn 100 #:thin 100))))<jupyter_output>cpu time: 300 real time: 300 gc time: 9<jupyter_text>Beyond uniform coin weightsSupposeyou know that the coin about to be flipped is weighted, so that it is more likely to land oneway than the other.  What you don't know is which side it favors, or to what degree. From your knowledge of the physics of flipping, though, extreme weights seem very unlikely.  To represent this state of knowledge, you need a distribution of *p(heads)* which concentrates its mass near 0.5.   Beta distributionsThere is a standard distibution family for just this case: the family of *beta* distributions.  A beta distribution assigns a probability density to every number in the interval $[0,1]$, and integrates to one over that interval.  A beta's shape is determined by two parameters, $\alpha$ and $\beta$, which act like "virtual observations".  That is, a preponderance of pseudo-observed ones (larger $\alpha$) moves the probability mass to the right, while a preponderance of pseudo-observed zeros (larger $\beta$) moves it to the left.  Here are a few beta distributions:<jupyter_code>(require "c3_helpers.rkt"
         racket/list
         math/distributions)<jupyter_output><empty_output><jupyter_text>$beta(1,1)$ is just a uniform distribution on $[0,1]$.<jupyter_code>(define (plot-beta alpha beta)
  (let ([xs (range 0.0001 0.9999 0.001)])
  (line-c3 xs 
           (map (lambda (x) (pdf (beta-dist alpha beta) x)) xs)
           #:xtickvalues '(0 1)
           #:showpoints #f)))

(plot-beta 1.5 2.5)

(pdf (beta-dist 1 1) 0.5)<jupyter_output><empty_output><jupyter_text>Like $beta(1,1)$, $beta(2,2)$ favors neither side of the interval, but, unlike $beta(1,1)$, it concentrates its probability mass near 0.5.  In general, if the $\alpha$ and $\beta$ parameters are given equal values, larger values will concentrate more sharply on 0.5.<jupyter_code>(show-beta 10 10)<jupyter_output><empty_output><jupyter_text>When $\alpha$ is big, relative to $\beta$, the probability mass shifts towards 1.<jupyter_code>(show-beta 4 1)<jupyter_output><empty_output><jupyter_text>Reversing the values of $\alpha$ and $\beta$ just flips the graph around the $x=0.5$ axis.$\alpha$ and $\beta$ values below 1 concentrate probability mass near the extremes:<jupyter_code>(show-beta 0.9 0.8)<jupyter_output><empty_output><jupyter_text>Assuming a moderate concentration of probability mass near 0.5, we can model 5 flips of the same weighted coin like so:<jupyter_code>(define weighted-flips
    (rejection-sampler
     
        ;; Generative model
        (define wt (beta 5 5))
        (deflazy result1 (bernoulli wt))
        (deflazy result2 (bernoulli wt))
        (deflazy result3 (bernoulli wt))
        (deflazy result4 (bernoulli wt))
        (deflazy result5 (bernoulli wt))
        (deflazy result6 (bernoulli wt))
     
        ;; Observations
        (observe result1 1)
        (observe result2 1)
        (observe result3 1)
        (observe result4 1)
        (observe result5 1)

        ;; Query expression
        (list wt (exact->inexact result6))))

   (sampler->mean weighted-flips 10000)<jupyter_output><empty_output><jupyter_text>Here we see that five 1s changes our best guess of the coin's weight to something like 0.67.You probably noticed the new deflazy [constructs] where the defines used to be.  These are necessary because observes combine with sampling expressions, under the hood, whereas defines are evaluated immediately. [...Ryan...]This happens to be one of those rare cases with a closed-form solution for the expected values of our query expressions.  The beta distribution is *conjugate* to the Bernoulli distribution, which means that the posterior distribution, given all of this observations, is also a beta distribution.  Remember where I said that $\alpha$ and $\beta$ behave like pseudo-observations?  Well, observations do too!  To get the posterior beta, we just add the number of 1s to alpha to define $\alpha'$, and the number of 0s to $\beta$ to define $\beta'$.  The posterior is then $beta(\alpha',\beta')$--here, $beta(10,5)$.A quick look at Wikipedia confirms that the mean of this posterior is $10/(10+5) = 2/3 \approx 0.66$, and the posterior predictive distribution of result6 is the same.[not budging enough after ten heads, leading to...][...more realistic, less mathematically convenient mixture distribution][...]Beyond BNs: indexed random variablesThere is a regular structure to the previous program that should be annoying to every programmer: all of those "define resultN" statements look the same.  All of the variables defined in those statements are drawn from the same distribution.  Can't we put them in an indexed set, such as an array or hash?  A popular Gamble idiom (borrowed from Church, but cleaned up) does something a tiny bit more abstract; we use defmem to draw a random *function* from objects of whatever type (integers, below) to draws from (bernoulli wt).  Conceptually, the entire (infinite) table of values is drawn at the time of definition, but the implementation is lazy; the random draw for a particular input to the random function is drawn the first time that input is seen, and memoized.  [What would be the most natural way to say "draw me a function at random from the set of functions that assign to each object in the range a draw from dist"?]<jupyter_code>;; Imports "range" and "for".
(require racket/list)

(define weighted-flips
   (rejection-sampler
   
    ;; Generative model
    (define wt (beta 5 5))
      (defmem (result flip-num)
        (bernoulli wt))
    
    ;; Observations
    (for ([x (range 1 6)])
          (observe (result x) 1))
    
    ;; Query expression
    (list wt (exact->inexact (result 6)))))

   (sampler->mean weighted-flips 10000)<jupyter_output><empty_output><jupyter_text>Metropolis-Hastings sampling also gets the right answer: (Hooray!) (define weighted-flips-mh   (mh-sampler       ;; Generative model    (define wt (beta 5 5))      (defmem (result flip-num)        (bernoulli wt))        ;; Observations    (for ([x (range 1 6)])          (observe (result x) 1))        ;; Query expression    (list wt (exact->inexact (result 6)))))   (sampler->mean weighted-flips-mh 1000 :burn 1000 :thin 50)  [We have to track this down and fix it: adaptive-drift-proposal: scale increased out of range]<jupyter_code>(expt 3 2)
(exact->inexact (expt (/ 51 49) 100))<jupyter_output><empty_output><jupyter_text>JunkSome inferences are easy...The "Asia" model is typical of an intuitive, causal style of modeling. It defines a joint probability distribution over risk factors, diseases, and symptoms (that is, it defines the probability of every distinct set of values for those variables) by breakingthe process in to small pieces, and it does so in a certain direction.  The pieces correspond to hypothesized "causal mechanisms", about which you can read more here (link).  For instance, disease probabilities depend directly on risk factors, because risk factors cause diseases.  Similarly, symptomprobabilities depend directly on diseases, because diseases cause symptoms.  Risk factors also causesymptoms, but only indirectly, via diseases.  A model of this form makes it easy to make one kind of conditional inference: forward, from causes to effects (note: the Bayes net formalism doesn't require a causal ordering, and arrows should not always beinterpreted as causal.  What makes inference easy, in any Bayes net, is query variables that are *desendants* of observed ones, whether those descendants are effects or not).Suppose we want to know the probability that Bob has tuberculosis, given that he has visited Asia.All we have to do is modify the model above to fix the values of $\verb+visit-to-asia+$ to $\verb+t+$. and more junkSo what?  You don't need the help ofprobability theory or computers to answer questions about a single coin flip, or even the long-run average of many coin flips.  If you have a little bit of experience with such things, you can probably even answer simple *conditional* questions like "what'sthe probability that the first of two flips was heads, given that at least one flip was heads.[more code...]WHERE?: In coin flipping, our ignorance of or inability to measure precisely enough the factors that determine the outcome is nearly complete, but for many other systems (e.g., the motion of aprojectile), we have pretty good deterministic models and measurements.  The probabilistic parts are there to captureminor deviations from our ideal model caused by such nuisances as a changeable atmosphere and imperfect sensors.  Random elements sweep under the rug all the factors that either aren't worth modeling, because they complicate the model without much improving predictions; or about which we can never know very much.     Intuition and case enumeration failWhen processes get just a little more complex than a pair of coin flips, though, unaided human minds either * fail to compute or intuit an answer:    * Q: If anywhere from one to six dice are selected, with equal probability, and the maximum across all dice, in one roll, is 5, what is the probability that the number of dice chosen is two?    * A: I need a computer for this.  Or at least lots of paper and time.* intuit wrong answers:     * Q: Are you having a hot streak at the roulette table?      * A: Yeah, Baby!  (wrong) Mathematical analysis failsSometimes we can use thecalculus of probability to answer questions our intution can't.  For instance, suppose .... Numerical approximations fail We can also use built-in programming language functions to predict the outcomes of more complicated processes,such as flipping a coin 100 times and outputing 1 if heads occurs at least 70 times.  What isthe probability of this event? If you said "about 0.0000016", give yourself a gold star (or give one to the software that told you the answer).Questions like that above (probability of 70+ heads in 100 fair coin flips) can be answered quickly, to a high degree of precision, by closed form math or fast, accurate numerical methods.  However, queries over more complex models have no such simple off-the-shelf solution methods. The progress of model development...[somewhere: Modeling the system under consideration--call it "s1"--is not the same thing as designing a system ("s2") to answer questions about s1.  s1 is usually what we care about, and s2 is what we have to muck with to get answers about s1. ][In these cases, writing down a model of the system you care about is only half (or a tenth of) the battle.  You also need to specify a method for getting an answer.]Simulation is useful when mathematical analysis runs out of steam.Inference: why we need probabilistic programmingMaking a model of how frequently diseases occur, and how often they produce various symptoms, is a fine exercise in formalization. First, write down a step-by-stepdescription of the process.  If it's too hard to do the math, or you are too lazy, then draw samples to get an estimate of how often it plays out this or that way.  Note, though,that, so far, we haven't done anything we couldn't have done in just about any existing programming language.  Why do we need anything new?  Mere simulation of a process is not enough, because **sampling from the possible executions of the process doesn'tsolve the kinds of problems people usually want to solve**.  Usually, we *observe* the values of some random variables in the model, and want to *infer* the probable values of some others.  For example, a doctorobserves a patient's symptoms and needs to infer the patient's probable disease status. The doctor *knows* that the process played out is some way such that particular symptomsoccurred, and wants to know something else about how it played out--what *unobservable* disease happenedto cause the symptoms.Some inferences are easy...The "Asia" model is typical of an intuitive, causal style of modeling. It defines a joint probability distribution over risk factors, diseases, and symptoms (that is, it defines the probability of every distinct set of values for those variables) by breakingthe process in to small pieces, and it does so in a certain direction.  The pieces correspond to hypothesized "causal mechanisms", about which you can read more here (link).  For instance, disease probabilities depend directly on risk factors, because risk factors cause diseases.  Similarly, symptomprobabilities depend directly on diseases, because diseases cause symptoms.  Risk factors also causesymptoms, but only indirectly, via diseases.  A model of this form makes it easy to make one kind of conditional inference: forward, from causes to effects (note: the Bayes net formalism doesn't require a causal ordering, and arrows should not always beinterpreted as causal.  What makes inference easy, in any Bayes net, is query variables that are *desendants* of observed ones, whether those descendants are effects or not).Suppose we want to know the probability that Bob has tuberculosis, given that he has visited Asia.All we have to do is modify the model above to fix the values of $\verb+visit-to-asia+$ to $\verb+t+$.Simple rejection sampling implements the meaning of conditioningThe simplest approximation is this: run the program many times, and throw out those runs that don't produce Bob's symptoms. The proportion of retained runs in which tuberculosis is $\verb+t+$ will be, on average,$p(tuberculosis|cough,x-ray)$.  Here we have * defined a program to   * draw a sample from the joint distribution, via a set of internal defines,  * define an output value ((map (lambda ...), based on that sample;* run that program 1000 times; and * displayed proportions of samples for which the three variables in the output expression take on the value t.  $\verb+asia+$ is a no-arg procedure (thunk) that draws a value for each random variable. sampler->mean is a convenience function that runs a sampling thunk a number of times and takes the element-wise means of the samples.Digression: visualizing resultsOne thing you will want to do quite often is *look at* probability distributions represented as set of samples.  Histograms are very useful for this.  For details, see the [plot](http://docs.racket-lang.org/plot/index.html?q=plot)library, but this examples gives an idea of the basic syntax: For some other small, special models, posterior distributionshave nice closed forms.  However, in the general case, these expressions involve very large sums and/orunsolvable integrals, and the only feasible approaches are approximate.  The next few sections discuss a few *Monte Carlo* approximations.  These are approximations that represent the posterior distributions of interest by sets of samples drawn from those distributions.  The first--rejection sampling--is a direct embodiment of the meaning of conditioning.  It is impractical for real problems, but useful as a conceptual aid, and as a check on the correctness of implementations of more sophisiticated algorithms.  The second is the (currently) most popular all purpose sampling method. Unfortunately, in many more complicated models nearly all runs will be rejected, because the observations are very improbable (in the case of continuous-valued observations, infinitelyso). So, while simple rejection sampling is an intuitive way to think about the meaning of conditionalinference, and sometimes a useful check on the correctness of implementations of other inference methods, itis rarely used in practice.  [Likelihood weighting]Metropolis-Hastings sampling is the Swiss army knife of inferenceWe solved the "Asia" diagnosis problem above using "enum".  We can also solve it approximately, using a general-purpose sampling algorithm:<jupyter_code>(require plot/pict
         racket/draw)
(plot (discrete-histogram `(#(dyspnea ,(list-ref asia-sample-means 0))
                         #(tuberculosis ,(list-ref asia-sample-means 1))
                         #(lung-cancer ,(list-ref asia-sample-means 2))))
      #:x-label "problems"
      #:y-label "probabilities of problems")
(require plot)
; A plotting convenience function
(define (show-beta alpha beta) 
  (plot (function 
         (lambda (x) (dist-pdf (beta-dist alpha beta) x)) 
         0.0001 
         0.9999 
         #:y-min 0
         #:y-max 4)))

(show-beta 1 1)
(define (asia-w-visit)
    (define visit-to-asia #t)
    (define smoker (flip 0.5))
    (define tuberculosis (if visit-to-asia (flip 0.05) (flip 0.01)))
    (define lung-cancer (if smoker (flip 0.1) (flip 0.01)))
    (define bronchitis (if smoker (flip 0.6) (flip 0.3)))
    (define tb-or-c (or tuberculosis lung-cancer))
    (define x-ray-abnormal (if tb-or-c (flip 0.98) (flip 0.05)))
    (define dyspnea (cond [(and tb-or-c bronchitis) (flip 0.9)]
                          [tb-or-c (flip 0.7)]
                          [bronchitis (flip 0.8)]
                          [else (flip 0.1)]))
    (map (lambda (x) (if x 1 0)) (list dyspnea tuberculosis lung-cancer)))

    (define asia-w-visit-sample-means (sampler->mean asia-w-visit 1000))
    
    (printf "dyspnea: ~a\ntuberculosis: ~a\nlung cancer: ~a\n" 
            (list-ref asia-w-visit-sample-means 0) 
            (list-ref asia-w-visit-sample-means 1) 
            (list-ref asia-w-visit-sample-means 2) )
(require plot/pict)
(plot (list (discrete-histogram `(#(dyspnea ,(list-ref asia-w-visit-sample-means 0))
                                  #(tuberculosis ,(list-ref asia-w-visit-sample-means 1))
                                  #(lung-cancer ,(list-ref asia-w-visit-sample-means 2)))
                                    #:skip 2.5 
                                    #:x-min 0
                                    #:label "Visit known")
            (discrete-histogram `(#(dyspnea ,(list-ref asia-sample-means 0))
                                  #(tuberculosis ,(list-ref asia-sample-means 1))
                                  #(lung-cancer ,(list-ref asia-sample-means 2)))
                                    #:skip 2.5 
                                    #:x-min 1
                                    #:label "Nothing known"
                                    #:color 2
                                    #:line-color 2))
          #:x-label "Problems"
          #:y-label "Problem probabilities")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06987069204728343, 'fraction_numerical': 0.009822964648426661, 'mean_word_length': 4.0488932474082375, 'pattern_counts': {'":': 1, '<': 94, '<?xml version=': 0, '>': 104, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '190417', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 9804, 'n_tokens_neox': 8960, 'n_words': 4876, 'path': 'sean_notebooks/gamble_intro/gamble_intro.ipynb', 'repo_name': 'ppaml-op3/iracket', 'size': '36068'}