Skip to content

Instantly share code, notes, and snippets.

@kkdd
Last active September 5, 2022 13:13
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 kkdd/99f4befb33eaff3d0aa31f99182cfd98 to your computer and use it in GitHub Desktop.
Save kkdd/99f4befb33eaff3d0aa31f99182cfd98 to your computer and use it in GitHub Desktop.
Drawing cumulative density function
<!DOCTYPE html>
<!-- ref: https://bl.ocks.org/jltran/bcd2a30fd9c08f8b9f87 -->
<meta charset="utf-8">
<style>
svg {
font: 10px sans-serif;
}
.axis path, .axis line {
fill: none;
stroke: #000;
shape-rendering: crispEdges;
}
</style>
<body>
<script src="https://d3js.org/d3.v7.min.js"></script>
<script>
//Set dimensions
const m = {top: 50, right: 50, bottom: 50, left: 50}
, h = 500 - m.top - m.bottom
, w = 960 - m.left - m.right
, numBins = 50
, n = 10000;
//Generate random numbers and calculate direct CDF
const norm = d3.range(n).map(d3.randomNormal(0, 1));
norm.sort((first, second) => first - second);
const cum = norm.slice(0, n - 1).map((w, i) => ({'x': w, 'y': (i+1)/n}));
//Calculate histogram and derived CDR
const xScale = d3.scaleLinear().domain([-3, 3]).range([0, w]);
const data = d3.histogram().domain(xScale.domain()).thresholds(xScale.ticks(numBins))(norm).filter(d => d.x0 < d.x1);
for(let i = 0, cum0 = 0; i < data.length; i++){
cum0 += data[i].length;
data[i]['cum'] = cum0 / n;
}
//Scales and axes
const yhist = d3.scaleLinear().range([h, 0]).domain([0, d3.max(data, d => d.length)]);
const ycum = d3.scaleLinear().domain([0, 1]).range([h, 0]);
const xAxis = d3.axisBottom(xScale);
const yAxis = d3.axisLeft(yhist);
const yAxis2 = d3.axisRight(ycum);
//Draw svg
const svg = d3.select("body").append("svg")
.attr("width", w + m.left + m.right)
.attr("height", h + m.top + m.bottom)
.append("g")
.attr("transform", "translate(" + m.left + "," + m.top + ")");
//Draw histogram
svg.selectAll("rect")
.data(data)
.enter()
.append("rect")
.attr("x", 1)
.attr("transform", d => "translate(" + xScale(d.x0) + "," + yhist(d.length) + ")")
.attr("width", d => xScale(d.x1) - xScale(d.x0) - 1)
.attr("height", d => h - yhist(d.length))
.style("fill", "#69b3a2")
.style("opacity", 0.2);
//Draw CDF's
svg.append("path")
.attr("d", d3.line().x(d => xScale(d.x1)).y(d => ycum(d.cum)).curve(d3.curveLinear)(data))
.attr("stroke", "#00ffff")
.attr("stroke-width", 7)
.attr("fill", "none");
svg.append("path")
.attr("d", (d3.line().x(d => xScale(d.x)).y(d => ycum(erfinv2(d.y))).curve(d3.curveLinear)(cum)))
.attr("stroke", "#ffc000")
.attr("stroke-width", 1.2)
.attr("fill", "none");
svg.append("path")
.attr("d", (d3.line().x(d => xScale(d.x)).y(d => ycum(d.y)).curve(d3.curveLinear)(cum)))
.attr("stroke", "#ff0000")
.attr("stroke-width", 1.2)
.attr("fill", "none");
//Draw axes
svg.append("g")
.attr("class", "x axis")
.attr("transform", "translate(0," + h + ")")
.call(xAxis);
svg.append("g")
.attr("class", "y axis")
.call(yAxis)
.append("text")
.attr("transform", "rotate(-90)")
.attr("y", 6)
.attr("dy", ".71em")
.style("text-anchor", "end")
.text("Count (Histogram)");
svg.append("g")
.attr("class", "y axis")
.attr("transform", "translate(" + [w, 0] + ")")
.call(yAxis2)
.append("text")
.attr("transform", "rotate(-90)")
.attr("y", 4)
.attr("dy", "-.71em")
.style("text-anchor", "end")
.text("CDF");
function erfinv2(x) {
return (erfinvapprox(x * 2 - 1) * 2 / Math.sqrt(Math.PI) + 1 ) / 2;
}
function erfinvapprox(x) {
// maximum relative error = .00013
const a = 0.147;
const log1 = Math.log(1 - x**2);
const b = 2 / (Math.PI * a) + log1/2
const sqrt1 = Math.sqrt(b**2 - log1/a)
return Math.sqrt(sqrt1 - b) * Math.sign(x)
}
//if (x <= 0) { return -999 } else if (x >= 1) {return 999}
</script>
</body>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment