function computeDensity(data, bins) {
  const [hist, binEdges] = computeHistogram(data, bins);
  // prevent dividing by 0
  const yScalingFactor =
    binEdges.length < 2 ? 1 : 1 / (binEdges[binEdges.length - 1] - binEdges[0]);
  return [hist.map(h => h * yScalingFactor), binEdges];
}
@param {Array} data - 1-D array of data points
@param {Number|Array} bins - number of bins, or bin edges
@return {Array<Array>} an array of 2 arrays:
0: density, histogram where each element is divided by the range of domain;
1: binEdges, an array of edge values (X values) of bins, length = hist.length + 1.