GCSR - Example 3.7

Attached is the Python source code for GCSR Example 3.7, the bioassay problem. The code is based on Helle Sorensen's code at KU. It samples from alpha, beta posteriors using grid approximation and displays a contour plot, and then, samples from the joint distribution and displays an X-Y plot of the sampled values. The plots duplicate Figure 3.4 (a) and (b) from GCSR.

Code

GCSR

I will publish sample code (Python), and solutions for selected questions found in Andrew Gelman's excellent book Bayesian Data Analysis 2nd Edition (GCSR). Each posting will be solution for one question, if required, will include mathematical formulas (using mathurl.com), and Python, R code. R code will be copy and paste from Gelman's own solutions PDF. We might also include code from Guilherme Rocha (another Bayesian Statistics lecturer), or from Jeff Hart (Texas A&M), or from Brian Junker (Carnegie Mellon University).

To install R on ubuntu, apt-get install r-base and r-base-dev.

You can run the R code from command line using "R -f [file]". If there are plots displayed in the R code, they will be dumped into a file called Rplots.pdf in your current directory.

Github

Sampling With Replacement Using Weights in Python

Here is the Python function corresponding to sample() call in R. We based it on the code here; only changed it so that the inputs use seperate weight and value vectors instead of one vector that has tuples of weight, value pairs.

import random

items = [(10, "low"),
       (100, "mid"),
       (890, "large")]

w = [10, 100, 890]
v = ["low", "mid", "large"]

def weighted_sample(ws, vs, n):
  total = float(sum(w for w in ws))
  i = 0
  w = ws[0]
  v = vs[0]
  while n:
      x = total * (1 - random.random() ** (1.0 / n))
      total -= x
      while x > w:
          x -= w
          i += 1
          w = ws[i]
          v = vs[i]
      w -= x
      yield v
      n -= 1

for i in weighted_sample(w, v, 500):
  print i

Bayes Election Prediction R Code, Gelman, Tokdar

Attached is the R code that calculates a prediction for 1992 presidential election based on methods found in Andrew Gelman's Bayesian Data Analysis book. The R code itself is written by Dr. Surya Tokdar; however the links at CMU that housed this and his other lecture related materials were all gone. I salvaged what I could from Google cache, removed all plotting, graphics related portions of the code, leaving only MCMC Gibbs calculations intact. Data was also missing at the CMU link, so I recreated both data files based on Andrew Gelman's original presidential.asc file. Still, the R chol() call (for Cholesky decomposition) gave errors, then I replaced all NA values with 0's. Code finally worked with this recent fix, and its results seem correct -- a Democrat win is predicted for most states at year 1992 which actually was the case. The code and data files can be found below. Let me know of your questions and comments.

Code