Skip to content

Instantly share code, notes, and snippets.

@EfratVil
Last active October 3, 2017 09:31
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 EfratVil/dedc46c4117ac11829e078673f076261 to your computer and use it in GitHub Desktop.
Save EfratVil/dedc46c4117ac11829e078673f076261 to your computer and use it in GitHub Desktop.
Mahalanobis Distance
<!DOCTYPE html>
<meta charset="utf-8">
<head>
<link href="https://fonts.googleapis.com/css?family=Raleway" rel="stylesheet">
<style>
body {
font-size: 20px;
font-family: 'Raleway', sans-serif;
}
</style>
</head>
<body>
<script src="https://d3js.org/d3.v4.min.js"></script>
<div style="margin-left: 0px; margin-top: 2px;">
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Select correlation level: <br/>
<span class="leftlabel">-1</span>
<input id ="range1" type="range" min="-99" max="99" value="70" style="width: 400px; margin-right: 10px;"/>
<span class="rightlabel">1</span>
<br />
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Current level: <span id="corr">0.7</span>
</div>
<script>
var margin = 25,
width = 400 - 2 * margin,
height = 400 - 2 * margin;
var svg = d3.select("body").append("svg")
.attr("id", "chart_svg")
.attr("width", width + 2 * margin)
.attr("height", height + 2 * margin)
.append("g")
.attr("transform", "translate(" + margin + " " + margin + ")");
var x = d3.scaleLinear()
.range([0, width]);
var y = d3.scaleLinear()
.range([height, 0]);
var xAxis = d3.axisBottom(x).ticks(12),
yAxis = d3.axisLeft(y).ticks(12 * height / width);
//---------------------------------- Create Bi-variate Normal -------------------------------------
// source: http://ryan-j-smith.github.io/D3-Examples/html/plotRandomCov.html
//-------------------------------------------------------------------------------------------------
var mu = [0, 0];
var sig = [[1, 0.7],
[0.7, 1]];
//2-D Cholesky decomposition of sigma
function chol2d(sig) {
var a, b, c;
a = Math.sqrt(sig[0][0]);
b = sig[0][1] / a;
c = Math.sqrt(sig[1][1] - b * b);
return [
[a, 0],
[b, c]
];
}
function randomData(samples) {
var data = [],
random = d3.randomNormal();
//Perform cholesky decomposition
var sqrtSig = chol2d(sig);
for (i = 0; i < samples; i++) {
var x = [random(), random()];
var y = [0, 0];
data.push({
x: mu[0] + sqrtSig[0][0] * x[0] + sqrtSig[0][1] * x[1],
y: mu[1] + sqrtSig[1][0] * x[0] + sqrtSig[1][1] * x[1]
});
}
return data;
}
function calculate_mahalanobis(ds) {
arr = [];
b = [];
for (var j = 0; j < ds.length; j++) {
b = [];
b.push(ds[j]["x"]);
b.push(ds[j]["y"]);
arr.push(b);
}
var columns = transpose(arr),
means = columns.map(mean),
invertedCovariance = invert(cov(columns, means));
var deltas = arr.map(function (row, i) {
return row.map(function (value, i) {
return value - means[i];
});
});
deltas.map(function (row, i) {
m = Math.sqrt(
multiply(multiply(row, invertedCovariance), row)
);
arr[i].push(m)
});
data = ds.map(function (d, i) {
return {
x: +d.x,
y: +d.y,
m: arr[i][2]
};
});
return data;
}
var ds = randomData(2000);
var data = calculate_mahalanobis(ds);
x.domain(d3.extent(data, function (d) { return d.x; })).nice();
y.domain(d3.extent(data, function (d) { return d.y; })).nice();
var maxDistance = 8;
var step = d3.scaleLinear()
.domain([1, 3])
.range([0,6]);
var color = d3.scaleLinear()
.domain([1, 1.5, 2, 2.5, 3, 3.5, 6])
.range(['#1a9850', '#66bd63', '#a6d96a', '#fdae61', '#f46d43', '#d73027', '#d73027'])
.interpolate(d3.interpolateHcl); //interpolateHsl interpolateHcl interpolateRgb
['#d73027', '#f46d43', '#fdae61', '#fee08b', '#ffffbf', '#d9ef8b', '#a6d96a', '#66bd63', '#1a9850']
svg.selectAll("circle")
.data(data)
.enter()
.append("circle")
.attr("class", "dot")
.attr("r", 3)
.attr("cx", function(d){return x(d.x);})
.attr("cy", function(d){return y(d.y);})
.attr("fill", function(d){
return color(d.m);
});
svg.append("g")
.attr("class", "x axis ")
.attr('id', "axis--x")
.attr("transform", "translate(0," + height + ")")
.call(xAxis);
svg.append("g")
.attr("class", "y axis")
.attr('id', "axis--y")
.call(yAxis);
//-------------------------------------------- math.js --------------------------------------------
// source: https://github.com/veltman/mahalanobis/blob/master/src/math.js
//-------------------------------------------------------------------------------------------------
function mean(arr) {
return sum(arr) / arr.length;
}
function sum(arr) {
return arr.reduce(function(a, b){
return a + b;
});
}
function isNumeric(n) {
return typeof n === "number" && !isNaN(n);
}
//------------------------------------------- matrix.js -------------------------------------------
// source: https://github.com/veltman/mahalanobis/blob/master/src/matrix.js
//-------------------------------------------------------------------------------------------------
function dot(a, b) {
if (a.length !== b.length) {
throw new TypeError("Vectors are of different sizes");
}
return sum(a.map(function(x, i){
return x * b[i];
}));
}
function multiply(a, b) {
var aSize = a.every(isNumeric) ? 1 : a.length,
bSize = b.every(isNumeric) ? 1 : b.length;
if (aSize === 1) {
if (bSize === 1) {
return dot(a, b);
}
return b.map(function(row){
return dot(a, row);
});
}
if (bSize === 1) {
return a.map(function(row){
return dot(row, b);
});
}
return a.map(function(x){
return transpose(b).map(function(y){
return dot(x, y);
});
});
}
function transpose(matrix) {
return matrix[0].map(function(d, i){
return matrix.map(function(row){
return row[i];
});
});
}
function cov(columns, means) {
return columns.map(function(c1, i){
return columns.map(function(c2, j){
var terms = c1.map(function(x, k){
return (x - means[i]) * (c2[k] - means[j]);
});
return sum(terms) / (c1.length - 1);
});
});
}
function invert(matrix) {
var size = matrix.length,
base,
swap,
augmented;
// Augment w/ identity matrix
augmented = matrix.map(function(row,i){
return row.slice(0).concat(row.slice(0).map(function(d,j){
return j === i ? 1 : 0;
}));
});
// Process each row
for (var r = 0; r < size; r++) {
base = augmented[r][r];
// Zero on diagonal, swap with a lower row
if (!base) {
for (var rr = r + 1; rr < size; rr++) {
if (augmented[rr][r]) {
// swap
swap = augmented[rr];
augmented[rr] = augmented[r];
augmented[r] = swap;
base = augmented[r][r];
break;
}
}
if (!base) {
throw new RangeError("Matrix not invertable.");
}
}
// 1 on the diagonal
for (var c = 0; c < size * 2; c++) {
augmented[r][c] = augmented[r][c] / base;
}
// Zeroes elsewhere
for (var q = 0; q < size; q++) {
if (q !== r) {
base = augmented[q][r];
for (var p = 0; p < size * 2; p++) {
augmented[q][p] -= base * augmented[r][p];
}
}
}
}
return augmented.map(function(row){
return row.slice(size);
});
}
//----------------------------------------- mahalanobis.js ----------------------------------------
// source: https://github.com/veltman/mahalanobis/blob/master/build/mahalanobis.js
//-------------------------------------------------------------------------------------------------
function mahalanobis(data, cov, loc_est) {
var deltas = row.map(function (d, i) {
return d - means[i];
});
}
d3.select("#range1").on("change", function () {
d3.select("#corr").html(d3.select("#range1").property("value") / 100);
sig = [[1, d3.select("#range1").property("value") / 100],
[d3.select("#range1").property("value") / 100, 1]];
ds = randomData(2000);
data = calculate_mahalanobis(ds);
x.domain(d3.extent(data, function (d) { return d.x; })).nice();
y.domain(d3.extent(data, function (d) { return d.y; })).nice();
var t = svg.transition().duration(750);
svg.select("#axis--x").transition(t).call(xAxis);
svg.select("#axis--y").transition(t).call(yAxis);
svg.selectAll(".dot").data(data)
svg.selectAll("circle").transition(t)
.attr("cx", function (d) { return x(d.x); })
.attr("cy", function (d) { return y(d.y); })
.attr("fill", function (d) {
return color(d.m);
});
});
</script>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment