Visualizing Bayes, Markov Chain Monte Carlo

Thanks to the excellent bayes.js library from Rasmus Bååth it's now possible to experiment with Bayesian statistics in JavaScript. We'll take advantage of that library in this series of posts, which demonstrate Bayesian statistics visually for anyone with a web browser.

This first post covers Markov Chain Monte Carlo (MCMC), algorithms which are fundamental to modern Bayesian analysis. MCMC is also critical to many machine learning applications. Since this is the first post, though, we'll start with a brief introduction to Bayesian statistics.

<!DOCTYPE html>
<html>
    <head>
        <meta charset="UTF-8">
        <title>bayes.js Demo</title>
        <script src="https://cdn.jsdelivr.net/gh/rasmusab/bayes.js/mcmc.js"></script>
        <script src="https://cdn.jsdelivr.net/gh/rasmusab/bayes.js/distributions.js"></script>
        <script src="https://d3js.org/d3.v3.min.js"></script>
        <link href='https://fonts.googleapis.com/css?family=Varela' 
              rel='stylesheet' type='text/css'>
        
        <style>
            body {
                color: #444;
                font-family: Varela,sans-serif;
            }
            .axis path, .axis line {
                fill: none;
                stroke: #888;
                shape-rendering: crispEdges;
            }
            .target {
                fill: #ca0000;
            }
        </style>
    </head>
    <body>
    
        <div id="graph"></div>
    
        <script>

// Standard D3.js scatterplot set up
            var margin = {top: 40, right: 40, bottom: 40, left: 40},
                width  = 636 - margin.left - margin.right,
                height = 476 - margin.top  - margin.bottom;

var svg = d3.select("#graph").append("svg")
                .attr("height", height + margin.left + margin.right)
                .attr("width", width + margin.top + margin.bottom);
            
            var chart = svg.append("g")
                .attr("transform",
                    "translate(" + margin.left + "," + margin.top + ")"
                );
                
            var xScale = d3.scale.linear()
                .range([0,width]);
            var yScale = d3.scale.linear()
                .range([height,0]);
            
            // Since we know that the actual mean and standard
            // devation are 100 and 10, respectively, we can set
            // the domains appropriately.
            xScale.domain([0, 200]);
            yScale.domain([0, 20]);
            
            var xAxis = d3.svg.axis()
                .scale(xScale)
                .orient("bottom");
            var yAxis = d3.svg.axis()
                .scale(yScale)
                .orient("left");
                
            chart.append("g")
                    .attr("class", "axis")
                    .attr("transform", "translate(0," + height + ")")
                    .call(xAxis)
                .append("text")
                    .attr("x", width)
                    .attr("y", -10)
                    .style("text-anchor", "end")
                    .text("μ");

chart.append("g")
                    .attr("class", "axis")
                    .call(yAxis)
                .append("text")
                    .attr("x", 15)
                    .attr("y", 6)
                    .style("text-anchor", "end")
                    .text("σ");

// Add a special target that the sampler should
            // converge around.
            chart.append("circle")
                    .attr("class", "target")
                    .attr("r", 6)
                    .attr("cx", xScale(100))
                    .attr("cy", yScale(10));
            
            chart.append("text")
                .attr("x", width)
                .attr("y", 20)
                .classed("notes", true)
                .style("text-anchor", "end")
                .text("Click to start");
            
            function addPoint(d, prev) {
                if (prev) { 
                    prev.transition().attr('fill-opacity', '0.1'); 
                }
                return chart.append("circle")
                    .attr('fill-opacity', '1')
                    .attr("r", 4)
                    .attr("cx", xScale(d.mu[0]))
                    .attr("cy", yScale(d.sigma[0]));
            }

// Set up the Bayesian analysis by defining the
            // parameters, the model, and the observations.
            // Create an MCMC sampler to evaluate the model.

// We want to model the observations as coming
            // from a Normal distribution with mean `mu`
            // and standard deviation `sigma`. So those
            // are our parameters. Both are real, and, of
            // course, we need to ensure that `sigma` is
            // positive.
            var params = {
                mu: {type: "real"},
                sigma: {type: "real", lower: 0}
            };
            
            // For the model itself, we first specify the
            // prior probabilities for the parameters. We
            // don't want to assume much about those parameters,
            // so we set the priors as follows:
            //
            // - `mu`: normal with mean of 0 and a huge
            //         standard deviation. Effectively,
            //         `mu` could be anything.
            // - `sigma`: uniform from 0 to 100. We're
            //            saying that `sigma` could be
            //            practically anything also, as
            //            long as it's positive.
            //
            // After specifying the prior probabilities,
            // we specify the likelihood. In this case,
            // we're simply modeling the observations as
            // data form a Normal distribution with a
            // a mean of `mu` (whatever that value turns
            // out to be) and a standard deviation of 
            // `sigma` (whatever it turns out to be).
            var model = function(state, obs) {
                var log_post = 0;
                // Priors
                log_post += ld.norm(state.mu, 0, 100);
                log_post += ld.unif(state.sigma, 0, 100);
                // Likelihood
                for(var i = 0; i < obs.length; i++) {
                    log_post += ld.norm(obs[i], state.mu, state.sigma);
                }
                return log_post;
            };
            
            // Here are the "observations." Don't tell
            // the sampler, but they were drawn from a 
            // normal distribution with mean 100 and
            // standard deviation of 10. We'll see how
            // close the sampler can get to those values.
            var obs = [
                120.6, 103.3, 108.1, 103.7, 96.35, 
                102.1, 96.23, 106.9, 111.2, 100.2
            ];
            
            var sampler =  new mcmc.AmwgSampler(params, model, obs);
            
            // Define the total number of samples
            // that we'll generate, as well as how
            // many of them we'll use as burn-in.            
            var numSamples = 1000, 
                burnin = 400;
            
            // Wait for a click to start the simulation
            d3.select('body').on('click', function() {

d3.select(".notes").text("");

Updated missing url https://rawgit.com/rasmusab/bayes.js/master/mcmc.js to https://cdn.jsdelivr.net/gh/rasmusab/bayes.js/mcmc.js

Updated missing url https://rawgit.com/rasmusab/bayes.js/master/distributions.js to https://cdn.jsdelivr.net/gh/rasmusab/bayes.js/distributions.js

https://rawgit.com/rasmusab/bayes.js/master/mcmc.js

https://rawgit.com/rasmusab/bayes.js/master/distributions.js

https://d3js.org/d3.v3.min.js