1. Inference in DAGs using rejection sampling

Here is a little Greek-gods-and-weather DAG written up in R, with an example of how to take “forward samples” from the full joint distribution on its variables.

flip = function(p) runif(1) < p
rain.sim = function(i) {
    flip = function(p) runif(1) < p
    Zeus.angry = flip(.3)
    Zeus.crying = flip(.2)
    thunder = (Zeus.angry && flip(.7)) || flip(.3)
    rain = (Zeus.angry && flip(.8)) || (Zeus.crying && flip(.99))
    forgetful.gardener = flip(.4)
    sprinkler = forgetful.gardener && flip(.7)
    wet.grass = (rain && flip(.8)) || (sprinkler && flip(.9)) || flip(.1)
    return(c(Zeus.angry = Zeus.angry, Zeus.crying = Zeus.crying, thunder = thunder, rain=rain,
             forgetful.gardener = forgetful.gardener, sprinkler = sprinkler, wet.grass = wet.grass))
}
sapply(1:10, FUN=rain.sim)
##                     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
## Zeus.angry         FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE
## Zeus.crying         TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE
## thunder             TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE
## rain                TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
## forgetful.gardener  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## sprinkler           TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE
## wet.grass          FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
##                    [,10]
## Zeus.angry         FALSE
## Zeus.crying        FALSE
## thunder            FALSE
## rain               FALSE
## forgetful.gardener  TRUE
## sprinkler           TRUE
## wet.grass          FALSE
  1. Work through the code and draw a DAG that encodes the relationships of immediate dependence and conditional independence between the variables. You can submit this part of the problem in any format you like (using fancy drawing software, pencil-and-paper, or anything in between - just make sure you put your name on the result!).
  2. What is the relationship between thunder and wet.grass? If you hear thunder while in bed at night, have you learned something about how likely it is that the grass will be wet in the morning? Why or why not?
  3. What if you also saw rain through the window - should learning about whether there was thunder still influence your best guess about the probability of wet grass?

Here’s an example of how to estimate probabilities from this DAG using Monte Carlo sampling:

query = function(variable.name, n.sims = 500) {
    sims = sapply(1:n.sims, FUN=rain.sim)
    proportion = length(which(sims[variable.name,] == TRUE))/n.sims
    return(proportion)
}
query('Zeus.angry')
## [1] 0.29
query('thunder')
## [1] 0.408
  1. For each variable, run the query() function to find out its approximate unconditional probability. (Or, more precisely, its approximate probability conditional on this model.) Make sure that you surround the variable name in quotes (single or double), both here and when working with the conditional.query() function below.

Most often, what we’re doing with a DAG is observing some ‘downstream’ variables (effects) and trying to infer the values of some ‘upstream’ variables (causes). This means that the unconditional probability is not really what’s of interest, but rather the conditional probability given our observations. Here’s an example of how to write a conditional query function that will estimate conditional probability using rejection sampling - a simple and reliable, but frequently inefficient, method of Monte Carlo approximation.

conditional.query = function(query.variable, conditioners, n.sims=1000) {
    samples = c()
    while (length(samples) < n.sims) {
      Zeus.angry = flip(.3)
      Zeus.crying = flip(.2)
      thunder = (Zeus.angry && flip(.7)) || flip(.3)
      rain = (Zeus.angry && flip(.8)) || (Zeus.crying && flip(.99))
      forgetful.gardener = flip(.4)
      sprinkler = forgetful.gardener && flip(.7)
      wet.grass = (rain && flip(.8)) || (sprinkler && flip(.9)) || flip(.1)
      truth.values = c(Zeus.angry = Zeus.angry, Zeus.crying = Zeus.crying, thunder = thunder, rain=rain,
                       forgetful = forgetful.gardener, sprinkler = sprinkler, wet.grass = wet.grass)
      query.index = which(query.variable == names(truth.values))
      conditioner.indices = which(names(truth.values) %in% conditioners)
      if (all(truth.values[conditioner.indices] == T)) { 
        samples = c(samples, truth.values[query.index])
      }
      # accept sample only if all conditioning statements are true
      # repeat until you have n.sims-many samples
    } 
    return(length(which(samples == T)) / n.sims) # proportion 'true' among accepted samples
}

  1. What is the estimated \(P(\mathbf{wet\ grass})\)? What about \(P(\mathbf{wet\ grass}|\mathbf{thunder})\)?, \(P(\mathbf{wet\ grass}|\mathbf{thunder},\mathbf{Zeus\ angry})\)? \(P(\mathbf{wet\ grass}|\mathbf{Zeus\ angry})\)? Explain, in terms of the conditional independence relationships encoded in the DAG you drew above, why the approximate probabilities are distributed in this way.

  2. How is the performance (runtime) of conditional.query() affected as you increase the number of items in the conditioners vector that you feed to this function? Using a for-loop (or something fancier if you want), test this question systematically for conditioners vectors with lengths from 0 to 7. Use system.time() to time the results.

  3. Explain, in terms of (a) the unconditional probability of the conditioning statements and (b) the implementation of conditional.query() using a while-loop, why the result in (1f) obtains.

Multivariate probability practice

  1. [Levy, exercise 3.1: Simpler formula for variance.] Show from the definition of the variance as \(\mathbf{Var}(X) ≡ \mathbb{E}[(X − \mathbf{E}(X))^2]\) that it can equivalently be written as \(\mathbf{Var}(X) = \mathbb{E}[X^2] − E[X]^2\), which we stated without proof in Section 2.9.2.

  2. [Levy, Exercise 3.2: Covariance of conditionally independent random variables.] Use linearity of the expectation [section 3.3.1] to prove that if two discrete random variables \(X\) and \(Y\) are conditionally independent given your state of knowledge, then \(\mathbf{Cov}(X, Y) = 0\). (Hint: you can rewrite \(P(X=x,Y=y) = \sum_{x,y} [x \times p(X = x) \times y \times p(Y = y)]\) as \([\sum_{x} x \times p(X = x)] \times [\sum_y y \times p(Y = y)]\), since \(X\) and \(p(X = x)\) are constant with respect to \(y\), and vice versa.)