Skip to content

Instantly share code, notes, and snippets.

@junosuarez
Created January 9, 2020 06:11
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 junosuarez/d1eec928f8840cac12aeda7cbe508dd6 to your computer and use it in GitHub Desktop.
Save junosuarez/d1eec928f8840cac12aeda7cbe508dd6 to your computer and use it in GitHub Desktop.
Approximating π by simulation
/*
Approximating π by simulation
see live on runkit at https://runkit.com/junosuarez/5e16c3e9ccfea2001ae2c8f5
In 1733 the French naturalist Georges Louis Leclerc, Comte de Buffon, posed and
solved the following problem in geometric probability: when a needle of length L is
dropped onto a wooden floor constructed from boards of width D (where D ≥ L),
what is the probability that it will lie across a seam between two boards?
Buffon determined that the probability is 2L/Dπ.
via https://www.maa.org/sites/default/files/s1000042.pdf
Can try to simulate this geometrically by trials of a random drop angle and origin and checking for
intersections of board lines
*/
function trial(L, D) {
// boards are parallel to the y axis. so it becomes a vertical line test.
// the needle is a line from A to B with length L
// assert precondition
if (!(D >= L)) {
throw new Error("D must be >= L");
}
// let our scale be +- 10 board widths
const scale = 10 * D;
// first we generate a random point A and a random slope m
const A = [Math.random(), Math.random()].map(n => n * scale);
const Ax = A[0];
const m = Math.random();
// then we solve the triangle to find B from L
// we know the hypotenuse and the slope, need to calculate opp and adj sides
// we can shortcut because we only need to know Bx, not B
const theta = 2 * Math.PI * m;
// by law of cosines, ratio between hypotenuse and adjacent side is cosine theta
const Δx = Math.cos(theta) * L;
const Bx = Ax + Δx;
// by Intermediate Value Theorem, we pass the vertical line test if Ax < Dn < Bx for any n in N
// console.log(Ax, "\t", Bx);
if (Ax === Bx) {
// if it's a vertical line, only pass if Ax = Dn for some n in N
return Ax % D === 0;
}
// naive iterate board widths, can be improved
// it doesn't matter which end of the line is A or B as long as the length is L
const min = Math.min(Ax, Bx);
const max = Math.max(Ax, Bx);
for (let x = 0; x < max; x += D) {
if (min < x && x < max) {
// console.log(min, "\t", x, "\t", max, "\t", true);
return true;
}
}
// console.log(min, "\t", "Dn", "\t", max, "\t", false);
return false;
}
function simulate(trials = 10) {
const L = 5;
const D = 10;
let pass = 0;
for (let i = 0; i < trials; i++) {
const result = trial(L, D);
if (result) {
pass++;
}
}
/*
Buffon's needle formula is p = 2L/Dπ
We simulated p, so we can solve for π:
π = 2L/Dp
*/
const p = pass / trials;
const π = (2 * L) / (D * p);
return { trials, p, π, e: Math.PI - π };
}
const results = [];
for (let j = 0; j < 3; j++) {
let result = simulate(1e5);
results.push(
Object.assign({ run: j, ["abs e"]: Math.abs(result.e) }, result)
);
}
// display a pretty table of results
const Table = require('easy-table')
const t = new Table();
results.forEach(r => {
Object.entries(r).forEach(([key, val]) => {
t.cell(key, val);
});
t.newRow();
});
t.sort(["abs e"]);
t.total("abs e", {
printer: function(val, width) {
return "Avg: " + val;
},
reduce: function(acc, val, idx, len) {
acc = acc + val;
return idx + 1 == len ? acc / len : acc;
}
});
console.log(t.toString());
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment