Animation visualizing the Bayesian Blocks dynamic programming algorithm for optimal histograms, inspired by Mike Bostock's Visualizing Algorithms. The algorithm is described in Studies in Astronomical Time Series Analysis. VI. Bayesian Block Representations.
For an example of Bayesian Blocks in action, see Bayesian Blocks Density Estimation.
| <!DOCTYPE html> | |
| <html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"> | |
| <head> | |
| <meta charset="utf-8" /> | |
| <title>Bayesian Blocks Algorithm</title> | |
| <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.6/d3.min.js" charset="utf-8"></script> | |
| <style> | |
| svg .histbar { | |
| stroke-width: 1; | |
| stroke: #000000; | |
| shape-rendering: crispEdges; | |
| } | |
| svg text { | |
| font: 16px sans-serif; | |
| } | |
| </style> | |
| </head> | |
| <body> | |
| <svg id="bbsvg" version="1.1" xmlns="http://www.w3.org/2000/svg" width="960" height="500" viewBox="0 0 960 500"> | |
| <g id="bars-g"></g> | |
| <g id="dots-g"></g> | |
| <g> | |
| <rect x="32" y="32" width="112" height="88" rx="16" ry="16" style="fill: #ffffcc; fill-opacity: 0.75;" /> | |
| <text id="outer-loop-text" x="64" y="64"></text> | |
| <text id="inner-loop-text" x="64" y="96"></text> | |
| </g> | |
| </svg> | |
| <script src="interactive-bayesian-blocks.js"></script> | |
| <script> | |
| "use strict"; | |
| // Generate data | |
| var nPoints = 100; | |
| var samples = [] | |
| for (var i=0; i<nPoints; ++i) { | |
| samples.push(d3.random.normal(64, 16)()); | |
| } | |
| // Set up plot | |
| var mysvg = document.getElementById("bbsvg"); | |
| var myBars = mysvg.getElementById("bars-g"); | |
| var myDots = mysvg.getElementById("dots-g"); | |
| var myILabel = mysvg.getElementById("outer-loop-text"); | |
| var myRLabel = mysvg.getElementById("inner-loop-text"); | |
| var mywidth = mysvg.viewBox.baseVal.width; | |
| var myheight = mysvg.viewBox.baseVal.height; | |
| var dotRadius = 8; | |
| var dotPadding = 4; | |
| var xScale = d3.scale.linear() | |
| .rangeRound([0, mywidth]) | |
| .domain([d3.min(samples) - 1, d3.max(samples) + 1]); | |
| var yScale = d3.scale.linear() | |
| .rangeRound([myheight-2*(dotRadius+dotPadding), 0]) | |
| .clamp(true) | |
| .domain([0, 0.8]); | |
| // Draw samples | |
| var dots = d3.select(myDots).selectAll("circle").data(samples); | |
| dots.enter().append("circle"); | |
| dots.attr("cx", function(d) { return xScale(d); }) | |
| .attr("cy", myheight - (dotRadius + dotPadding)) | |
| .attr("r", dotRadius); | |
| // Construct algorithm | |
| var ncp_prior = 0.1; | |
| var alg = new BayesianBlocksAlgorithm(samples, ncp_prior); | |
| // Assign colors to ops | |
| var colors = []; | |
| colors[alg.states.FIT_NULL] = "#cccccc"; | |
| colors[alg.states.FIT_INIT] = "#cccc33"; | |
| colors[alg.states.FIT_IMPROVED] = "#33cc33"; | |
| colors[alg.states.FIT_REGRESSED] = "#cc3333"; | |
| colors[alg.states.FIT_SELECTED] = "#3333cc"; | |
| function updateViz() { | |
| // Update histogram bins | |
| // Use bin end as key for object constancy | |
| var bars = d3.select(myBars).selectAll("rect").data(alg.bins, function(d) { return d.end; }); | |
| bars.enter().append("rect") | |
| .classed("histbar", true) | |
| .attr("x", function(d) { return xScale(d.start); }) | |
| .attr("width", function(d) { return 0; }) | |
| .attr("y", function(d) { return yScale(d.height); }) | |
| .attr("height", function(d) { return yScale(0) - yScale(d.height); }) | |
| .style("fill", colors[alg.states.FIT_NULL]); | |
| bars.transition() | |
| .attr("x", function(d) { return xScale(d.start); }) | |
| .attr("width", function(d) { return xScale(d.end) - xScale(d.start); }) | |
| .attr("y", function(d) { return yScale(d.height); }) | |
| .attr("height", function(d) { return yScale(0) - yScale(d.height); }) | |
| .transition() | |
| .style("fill", colors[alg.lastOp]) | |
| .transition() | |
| .style("fill", colors[alg.states.FIT_NULL]); | |
| bars.exit().transition() | |
| .attr("width", 0) | |
| .remove(); | |
| // Update status text | |
| d3.select(myILabel).text("R = " + alg.lastI); | |
| d3.select(myRLabel).text("r = " + ((alg.lastR >= 0) ? alg.lastR : alg.bestR)); | |
| } | |
| // Animate the algorithm | |
| (function animateAlg() { | |
| if (!alg.isFinished()) { | |
| alg.step(); | |
| updateViz(); | |
| setTimeout(animateAlg, 1000); | |
| } | |
| })(); | |
| </script> | |
| </body> | |
| </html> |
| "use strict"; | |
| var BayesianBlocksAlgorithm = function(x, ncp_prior) { | |
| // Sort the input data | |
| var samples = x.slice().sort(d3.ascending); | |
| // Compute Voronoi cells for data | |
| var boundaries = [samples[0]]; | |
| for (var i=1; i<samples.length; ++i) { | |
| boundaries.push(0.5*(samples[i] + samples[i-1])); | |
| } | |
| boundaries.push(samples[samples.length-1]); | |
| /* Public output/state of the algorithm */ | |
| this.states = { | |
| FIT_NULL: 0, | |
| FIT_INIT: 1, | |
| FIT_IMPROVED: 2, | |
| FIT_REGRESSED: 3, | |
| FIT_SELECTED: 4 | |
| }; | |
| this.bins = []; | |
| this.lastOp = this.states.FIT_NULL; | |
| this.lastR = -1; | |
| this.bestR = NaN; | |
| this.lastI = 0; | |
| /* Private state of the algorithm */ | |
| var fits = [] | |
| var lasts = [] | |
| var bestFit = NaN; | |
| /* Private methods */ | |
| function blockFitness(r, i) { | |
| var n = i + 1 - r; | |
| var t = boundaries[i + 1] - boundaries[r]; | |
| return n*(Math.log(n) - Math.log(t)); | |
| } | |
| function mkBins(i, r) { | |
| function blockHeight(r, i) { | |
| var n = i + 1 - r; | |
| var t = boundaries[i + 1] - boundaries[r]; | |
| return n/t; | |
| } | |
| var left = r; | |
| var right = i+1; | |
| var ans = []; | |
| while (right > 0) { | |
| ans.push({ | |
| "start": boundaries[left], | |
| "end": boundaries[right], | |
| "height": blockHeight(left, right-1) | |
| }); | |
| right = left; | |
| left = (left > 0) ? lasts[left-1] : 0; | |
| } | |
| return ans.reverse(); | |
| } | |
| /* Privileged methods */ | |
| this.isFinished = function() { | |
| return fits.length >= samples.length; | |
| }; | |
| this.step = function() { | |
| var i = fits.length; | |
| if (this.lastR == 0) { | |
| // We have finished searching for the best r; select and show our choice | |
| fits.push(bestFit); | |
| lasts.push(this.bestR); | |
| this.lastR = -1; | |
| this.bins = mkBins(i, lasts[i]); | |
| this.lastOp = this.states.FIT_SELECTED; | |
| } else if (this.lastR == -1) { | |
| // Add the next data cell and initialize the search for r | |
| var r = i; | |
| bestFit = ((r > 0) ? fits[r-1] : 0.0) + blockFitness(r, i) - ncp_prior; | |
| this.bestR = r; | |
| this.bins = mkBins(i, r); | |
| this.lastOp = this.states.FIT_INIT; | |
| this.lastR = r; | |
| } else { | |
| // Evaluate next r | |
| var r = this.lastR - 1; | |
| var fit = ((r > 0) ? fits[r-1] : 0.0) + blockFitness(r, i) - ncp_prior; | |
| this.bins = mkBins(i, r); | |
| if (fit > bestFit) { | |
| bestFit = fit; | |
| this.bestR = r; | |
| this.lastOp = this.states.FIT_IMPROVED; | |
| } else { | |
| this.lastOp = this.states.FIT_REGRESSED; | |
| } | |
| this.lastR = r; | |
| } | |
| this.lastI = i; | |
| }; | |
| }; |
