Skip to content

Instantly share code, notes, and snippets.

@eliangcs
Last active January 18, 2016 06:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save eliangcs/6e8b45f88fd3767363e7 to your computer and use it in GitHub Desktop.
Save eliangcs/6e8b45f88fd3767363e7 to your computer and use it in GitHub Desktop.
Integrating sine with Monte Carlo / Metropolis algorithm
<!DOCTYPE html>
<meta charset="utf-8">
<style>
.svg {
border: 1px solid #000;
}
.result {
color: #05c;
}
</style>
<body>
<div class="result"></div>
<script src="//d3js.org/d3.v3.min.js"></script>
<script>
var width = 540,
height = 320,
xmin = 0,
xmax = Math.PI,
ymin = 0,
ymax = 1.0,
xscale = d3.scale.linear(),
yscale = d3.scale.linear(),
sine = d3.svg.line(),
svg = d3.select('body').append('svg'),
decor = svg.append('g'),
graph = svg.append('g'),
sampleSize = 100,
sampleCount = 0,
delay = 10,
integral = 0.0;
function reset() {
integral = 0.0;
sampleCount = 0;
graph.selectAll('*').remove();
integralResult.text('-');
}
xscale.domain([xmin, xmax])
.range([0, width]);
yscale.domain([ymin, ymax])
.range([height, 0]);
svg.attr({
width: width,
height: height,
class: 'svg'
});
function draw(x, y, area) {
var h = y;
var w = area / h;
x = xscale(x - w * 0.5);
y = yscale(y);
w = xscale(w);
graph.append('rect').attr({
x: x,
y: y,
width: w,
height: height - y,
fill: '#06c',
opacity: 0.2
});
}
// The function f(x) we're integrating
function f(x) {
if (x >= 0.0 && x <= Math.PI) {
return Math.sin(x);
}
return 0.0;
}
// The probability density of a sample x
function p(x) {
// Uniform distribution over [0, PI]
return 1.0 / Math.PI;
}
// Given a current location x, propose a new location x2
function mutate(x) {
// Uniform sampling over [0, PI]
return Math.random() * Math.PI;
}
// Given the current location x, what's the probability that the algorithm
// accepts the new location x2
function accept(x, x2) {
// return 1;
return Math.min(1.0, p(x2) / p(x));
}
// Metropolis algorithm
function sample(x) {
var x2 = mutate(x);
var a = accept(x, x2);
if (Math.random() < a) {
x = x2;
}
var y = f(x);
var d = y / (p(x) * sampleSize);
integral += d;
draw(x, y, d);
sampleCount++;
if (sampleCount < sampleSize) {
setTimeout(function() {
sample(x);
}, delay);
} else {
d3.select('.result').text('Integral: ' + integral);
}
}
var x0 = Math.random() * Math.PI;
sample(x0);
</script>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment