The avatar for @arnavm
arnavm
Published unlisted
Edited
Dec 7, 2021
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
function p0_prior(p0, length) {
// Empirical prior, parametrized by the false alarm probability ``p0``
// See eq. 21 in Scargle (2012)
return 4 - Math.log(73.53 * p0 * Math.pow(length, -0.478))
}
Insert cell
function objective_OP(T_k, N_k) {
return math.dotMultiply(N_k, math.subtract(math.log(N_k), math.log(T_k)));
}
Insert cell
function objective_PELT(T_k, N_k) {
return math.multiply(-1, objective_OP(T_k, N_k));
}
Insert cell
function* bayesian_blocks(data, mode, p0, stepby) {
// Sort input array
const sorted_data = data.slice().sort((a, b) => a - b);
const N = sorted_data.length;
// Calculate edges of each cell
const edges = getEdges(sorted_data)
// Calculate length of each block
const block_length = math.subtract(sorted_data[N - 1], edges);
// Arrays to store best configurations
var best = math.zeros(N)._data;
var last = math.zeros(N)._data;
// Array to save changepoints
var change_points = math.zeros(N)._data;
// Set mode
if (mode === undefined) mode = "OP";
// Get prior
if (p0 === undefined) p0 = 0.05;
var prior = p0_prior(p0, N);
// Get counts of events in individual cells
// This holds (virtually) true here since
// x coordinates are random reals.
const x = math.ones(N)._data;
// Additional bookkeeping for PELT
if (mode == "PELT") {
var i_min;
var t_max;
var T_k_subset;
var N_k_subset;
var peltFitness;
var candidates;
var pruned;
var K = 0;
var unpruned = [];
var lastPruned = -1;
best = math.zeros(N + 1)._data;
}
// Variables needed for iteration
var T_k, N_k_reversed, N_k, fit_vec, A_R, i_max, calculations = 0;
// Iterate over the cells
for (let i = 0; i < N; i++) {
if (mode == "OP") {
// Width of each block
T_k = math.subtract(block_length.slice(0, i + 1), block_length[i + 1]);
// Number of elements in each block
N_k_reversed = cumulativeSum(x.slice(0, i + 1).reverse());
N_k = N_k_reversed.reverse();
// Calculate objective function
fit_vec = objective_OP(T_k, N_k);
// A_R is the fitness values adjusted by the prior
A_R = math.subtract(fit_vec, prior);
A_R = Array.prototype.concat(
A_R.slice(0, 1),
math.add(A_R.slice(1, N), best.slice(0, i))
);
// Find index where fitness is maximized
i_max = argMax(A_R);
// Save the location of the last block
// and the best fitness found so far
last[i] = i_max;
best[i] = A_R[i_max];
if (stepby == "Cells") yield last[i];
}
if (mode == "PELT") {
unpruned = math.concat(unpruned, [i]);
T_k_subset = math.subset(block_length, math.index(unpruned));
if (unpruned.length == 1) T_k_subset = [T_k_subset];
T_k = math.subtract(T_k_subset, block_length[i + 1]);
N_k_subset = math.subset(x, math.index(unpruned));
if (unpruned.length == 1) N_k_subset = [N_k_subset];
N_k_reversed = cumulativeSum(N_k_subset.reverse());
N_k = N_k_reversed.reverse();
fit_vec = objective_PELT(T_k, N_k);
A_R = math.add(fit_vec, prior);
A_R = math.add(A_R, best.slice(lastPruned + 1, i + 1));
// Find index where cost is minimized
i_min = argMin(A_R);
last[i] = unpruned[i_min];
best[i + 1] = A_R[i_min];
peltFitness = math.add(best.slice(lastPruned + 1, i + 1), fit_vec);
candidates = indexWhere(N_k.slice(0, N_k.length - 1), (x) => x >= 2)
if (candidates.length > 0) {
pruned = indexWhere(math.larger(math.add(peltFitness, K), best[i + 1]), (x) => x == true);
if (pruned.length > 0) {
i_max = math.max(pruned);
if (arraysEqual(pruned, math.range(0, i_max + 1)._data)) i_max = math.max(pruned);
else i_max = -1;
if (i_max >= 0) {
t_max = unpruned[i_max];
unpruned = unpruned.slice(i_max + 1);
lastPruned = t_max;
}
}
}
if (stepby == "Cells") yield [last[i], lastPruned];
}
if (stepby == "Points") yield {"x": sorted_data[i], "delay": fit_vec.length};
}
// Find changepoints by iteratively peeling off the last block
var i_cp = N;
var ind = N;
while (true) {
i_cp--;
change_points[i_cp] = ind;
if (ind == 0) break;
ind = last[ind - 1];
}
change_points = change_points.slice(i_cp, N);
yield edges.get(change_points);
}
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
Insert cell
d3 = require("d3@5")
Insert cell
import {dodger} from "@d3/d3-random"
Insert cell
import {button, radio, slider} from "@jashkenas/inputs"
Insert cell
math = require('mathjs@6.0.4/dist/math.min.js')
Insert cell
import {dc} from "@gjmcn/data-cube-array-oriented-javascript"
Insert cell
import {reference, data} from "@mbostock/reference"
Insert cell

One platform to build and deploy the best data apps

Experiment and prototype by building visualizations in live JavaScript notebooks. Collaborate with your team and decide which concepts to build out.
Use Observable Framework to build data apps locally. Use data loaders to build in any language or library, including Python, SQL, and R.
Seamlessly deploy to Observable. Test before you ship, use automatic deploy-on-commit, and ensure your projects are always up-to-date.
Learn more