Introduction

• PPAML Summer School 2016.
• Language intro (DIPPL).
• Language intro (AgentModels).
• Bayesian Conitive Modeling, A Practical Course; M. D. Lee and E-J. Wagenmakers, Cambridge University Press, 2013. (BCM Book)

Simple Start (Sum)

Goal: If U1 and U2 are two independent variables with uniform distribution over the integers 1, 2, and 3, what is the distribution of U1 + U2 ?

var main = function () {

  var u1 = uniformDraw([1,2,3])
  var u2 = uniformDraw([1,2,3])
  
  return u1 + u2
}

var dist = Infer(
  {method: 'enumerate'},
  main);


print("Pr(x=1) = " + Math.exp(dist.score(1)).toFixed(4))
print("Pr(x=2) = " + Math.exp(dist.score(2)).toFixed(4))
print("Pr(x=3) = " + Math.exp(dist.score(3)).toFixed(4))
print("Pr(x=4) = " + Math.exp(dist.score(4)).toFixed(4))
print("Pr(x=5) = " + Math.exp(dist.score(5)).toFixed(4))
print("Pr(x=6) = " + Math.exp(dist.score(6)).toFixed(4))

viz.table(dist);
viz.auto(dist);

display("Most likely value and its probability: ")
print(dist.MAP().val.toFixed(4) + " its Pr = " + (Math.exp(dist.MAP().score).toFixed(4)))


var variance = function(a) {
  var m = expectation(a);
  return expectation(a, function(x) { return Math.pow(x - m, 2) })
}

print("Mean: " + expectation(dist).toFixed(4))
display("Variance: " + variance(dist).toFixed(4))

Consider: What are the distributions of U1 - U2, U1 * U2, and U1 / U2? What changes if the distribution of U2 is uniformDraw([0,1,2])?

Consider: The distribution returns the log-probability of elements using dist.score(val). To get the probability, take the exponent. To print, convert to a fixed point value (here with four decimals).

Simple Start (Comparison)

Goal: If U1 and U2 are two independent variables with uniform distribution over the integers 1, 2, and 3, what is the distribution of U1 == U2 ?

var main = function () {

  var u1 = uniformDraw([1,2,3])
  var u2 = uniformDraw([1,2,3])
  
  return u1 == u2
}

var dist = Infer(
  {method: 'enumerate'},
  main);

viz.table(dist);
//viz.auto(dist);

Consider: How can you simplify the above program to draw only a single sample (from another distribution)?

Consider: How does the inference results change if we instead use this inference method: var dist = Infer( {method: 'MCMC', samples: 100}, main); (try running several times).

Consider: How much do the run-to-run differences change if we increase the number of samples to 1000 or 10000?

Several More Simple Models (AgentModels)

WebPPL is a functional subset of JavaScript, with several probabilistic constructs. To get accustomed with the syntax, read and try out the examples from this page.

Consider:

• How to filter only numerical array elements? (e.g., this is useful for data cleaning)
• What do the functions map, filter, and zip do? (How do they relate to loops?) See all built-in functions here.
• How is recursion used to implement a geometric distribution?
• See different ways to plot data in 1D or 2D.

Optionally, see also this introduction.

Another Simple Model

Goal: We have three independent coin flips, assuming they are fair. We observed at least two heads as the outcome. What is the (posterior) probability that the first coin flip resulted in a head? What is the posterior expected value of that flip?

var model = function() {
  var c1 = flip(0.5);
  var c2 = flip(0.5);
  var c3 = flip(0.5);

  condition (c1+c2+c3 >= 2)
  
  return c1;
}

var dist = Infer( {method: 'enumerate'}, model);

viz.table(dist)

print("Expected value: " + expectation(dist).toFixed(3))


viz.auto(Infer({method: 'MCMC', samples: 10000}, model));

Consider:How is the expectation of the Bernoulli variable computed? What is it equal to?

Consider:How much does the prediction change when we change the prior distribution of c1 to 0.6?

Consider:Write a function that prints the expectation for every change of the prior probability of c1 from 0.4 to 0.6. Hint: use the function 'map'.

Simple Rate Inference (Beta-Bernoulli)

Goal: Beta distribution is a common prior for the Bernoulli distribution -- because its support is [0,1], i.e., the range of the parameter of the Bernoulli distribution. We want to infer the underlying probability of the coin that is tossed three times:

  var main = function () {
  var theta = beta(1,1)
  
  var k1 = flip(theta)
  var k2 = flip(theta)
  var k3 = flip(theta)

  condition (k1+k2+k3 == 2)  
  
  return theta
}

var dist = Infer(
  {method: 'MCMC', samples: 100000},
  main);
viz.auto(dist);

Consider: How would you modify the example to model tossing 5 coins? 7 coins?

Consider: Why we cannot use the inference method 'enumerate' for this example?

Consider: Try plotting the prior (beta(1,1)). You can do it by using the inference method 'forward'. What distribution does it look like? (see also discussion later about prior redictive and posterior redictive).

Simple Rate Inference (Beta-Binomial)

Goal: Infer the underlying success rate for a binary process (a binomial distribution represents a sum of independent coin tosses). We want to determine the posterior distribution of the rate theta.

  var main = function () {
  var n = 10 //100
  var kk = 5 //50
  var theta = beta(1,1)
  
  var k = binomial(theta, n)
  
  condition (k == kk)  
  
  return theta
}

var dist = Infer(
  {method: 'MCMC', samples: 10000},
  main);
viz.auto(dist);

var main = function () {
  var n1 = 10
  var n2 = 10
  var theta1 = beta(1,1)
  var theta2 = beta(1,1)

  
  var k1 = binomial(theta1, n1)
  var k2 = binomial(theta2, n2)
  
  condition (k1 == 5)
  condition (k2 == 7)
  
  var delta = theta1 - theta2
  
  return delta
}

var dist = Infer(
  {method: 'MCMC', samples: 100000},
  main);

viz.auto(dist);

// same as before
var main = function () {
  var n1 = 10
  var n2 = 10
  var theta1 = beta(1,1)
  var theta2 = beta(1,1)

  
  var k1 = binomial(theta1, n1)
  var k2 = binomial(theta2, n2)
  
  condition (k1 == 5)
  condition (k2 == 7)
  
  var delta = theta1 - theta2
  
  return [theta1, theta2] // note: here we return the joint distribution
}

var dist1 = Infer(
  {method: 'MCMC', samples: 10000},
  main);

viz.auto(dist1);
viz.marginals(dist1); // we can plot the marginals directly

  var main = function () {
  var n = 10 
  
  // prior: 
  var theta = beta(1,1)
  
  // posterior: 
  var k = binomial(theta, n)
  
  condition (k == 5)  

  // prior predictive: 
  var thetaprior = beta(1,1)
  var priorpredictive = binomial(thetaprior, n)

  // posterior predictive:
  var postpredictive = binomial(theta, n)

  return [theta, priorpredictive, postpredictive]

}

var dist = Infer(
  {method: 'MCMC', samples: 10000, burn: 1000, lag: 2},
  main);
viz.marginals(dist);

var trueObservations = [false, true, false];

// The state transition (unobserved):
// transition to the other state with the probability 0.3
var transprob = 0.3
var transition = function(s) {
  return s ? flip(1-transprob) : flip(transprob)
}

// Noisy observation:
var noiselvl = 0.1
var observe = function(s) { 
  var changeobservation = flip(noiselvl)
  return changeobservation ? !s : s 
}

var hmm = function(n) {
  var prev = (n > 1) ? // note recursion !!
    hmm(n - 1) :       // <---- here!
  	{states: [false], observations: []} // initial state

  // moving to the next state (the system itself is probabilistic):
  var newState = transition(prev.states[prev.states.length - 1]);

  // observe the noisy data
  var newObs = observe(newState);

  // return the trace of the states and the oboservations (we compute the probability of those) 
  return {
    states: prev.states.concat([newState]),
    observations: prev.observations.concat([newObs])
  };
}

// run hmm with the three observations
var hmmrunner =   function() {
    var r = hmm(trueObservations.length);
    map2(function(o,d) {condition(o===d)}, r.observations, trueObservations)
    return r.states;
}

var stateDist = Infer({method: 'enumerate'}, hmmrunner)

viz.table(stateDist)