Public
Edited
Jul 24, 2020
1 star
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
function randomVonMises(mu, kappa) {
if ((kappa = +kappa) <= 0) {
return () => (1 - 2 * Math.random()) * Math.PI;
}
const s = 0.5 / kappa,
r = s + Math.hypot(1, s);
return function() {
do {
var z = Math.cos(Math.PI * Math.random()),
d = z / (r + z),
u = Math.random();
} while (1 - d * d <= u && (1 - d) * Math.exp(d) < u);
var f = (1 + r * z) / (r + z);
return wrap(mu + Math.sign(Math.random() - 0.5) * Math.acos(f));
}
}
Insert cell
Insert cell
Insert cell
Insert cell
function randomVonMisesFisher(n, kappa) {
const m = n - 1,
h = Math.hypot(2 * kappa, m),
b = m / (h + 2 * kappa),
a = (2 * kappa + h + m) / 4,
d = m * (1 - Math.log(m)),
sb = symmetricBeta(m / 2),
Z = d3.randomNormal();
return function() {
do {
var y = sb(),
denom = 1 - (1 - b) * y,
t = 2 * a * b / denom;
} while (m * Math.log(t) - t + d < Math.log(Math.random()));
do {
var vec = d3.range(m).map(Z),
ssq = vec.reduce((a, c) => a + c * c, 0);
} while (ssq === 0);
var w = (1 - (1 + b) * y) / denom,
mult = Math.sqrt((1 - w * w) / ssq);
vec = vec.map(x => x * mult);
vec.push(w);
return vec;
}
}
Insert cell
// Ulrich also provides the following generator for symmetric beta variates
// that is faster than blindly computing two gamma variates.
function symmetricBeta(a) {
if ((a = +a) === 0.5) {
return () => (Math.cos(Math.PI * Math.random()) + 1) / 2;
}
if (a === 1) {
return () => Math.random();
}
const b = 1 / (a - 0.5);
return function() {
do {
var u = 2 * Math.random() - 1,
v = 1 - Math.random(),
s = u * u + v * v;
} while (s > 1);
return u * v * Math.sqrt(1 - s ** b) / s + 0.5;
}
}
Insert cell
Insert cell
function mleKappa(r, p = 2) {
var k = r * (p - r * r) / (1 - r * r),
ak, delta = 1;
while (Math.abs(delta) > 1e-12) {
ak = besselRatio(p / 2 - 1, k);
delta = (ak - r) / (1 - ak * ak - (p - 1) / k * ak);
k -= delta;
}
return k;
}
Insert cell
// Computes I_{nu+1}(x) / I_{nu}(x).
// Amos (1974), Computation of modified Bessel functions and
// their ratios, Mathematics of Computation, 28 (125), 235-251
function besselRatio(nu, x) {
var l = Array(),
res = 2;
for (let k = 0;; k++) {
let s = nu + k + 0.5;
l.push(x / (s + Math.hypot(s + 1, x)));
for (let m = k; m > 0; m--) {
let r = Math.sqrt(l[m] / l[m-1]),
s = nu + m;
l[m-1] = x / (s + Math.hypot(s, x * r));
}
if (Math.abs(l[0] - res) < 1e-16) {
return res;
}
res = l[0];
}
}
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell

Purpose-built for displays of data

Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
Learn more