Added density reducer to bin/group/hexbin

[wirhabenzeit]

This is a pull request addressing #1940, see also the discussion in https://talk.observablehq.com/t/normalised-histogram-in-observable-plot/8576

Basically I added a new reduceDensity reducer which computes the density per series/facet. This is slightly different from the reduceProportion (with or without the facet scope) reducer (when not supplied with a value)

• reduceDensity does normalise by the bin size
• reduceDensity normalises by group rather than globally or by facet.

I forked the notebook https://observablehq.com/@fil/plot-normalized-histograms here https://observablehq.com/d/fb0d876105777d59 to illustrate the functionality. I also added a test plot in test/plots/density-reducer.ts

The main change to existing code is that in src/transforms/group.js, src/transforms/bin.js and src/transforms/hexbin.js need to call the reducer scope on a per group level as in

for (const o of outputs) o.scope("group", I);

I noticed that this new density reducer can also be used as a replacement for proportion-facet in

• test/plots/athletes-sport-weight.ts
• test/plots/hexbin-r.ts

where (at least to me) it also makes sense semantically.

Among the tests this would leave

• test/plots/penguin-species-island-relative.ts

as the only place where proportion-facet is used. Here density does not makes sense on a group level due to the grouping by fill.

One could also consider a more customisable reducer allowing to specify the normalisation scope, but I could not come up with a satisfying syntax. Something like

..., y: {reducer: "density", scope: "group"}, ...

feels a bit too verbose given that most reducers do not have/need any scope.

[Fil]

This is looking great!

Suggestions:

I think we can extend the concept to weighted density when the original channel is a value?

The documentation needs to clarify what the reducer does when you're not reducing y in binX or x in binY.

If we're binning in two dimensions, and using this reducer to inform r, then the "area" is defined by x and y and is not 1, and the results are comparable to what happens with hexbin.

The documentation should also clarify that this normalization is applied on each series (all values with the same z in the same facet).

[Fil]

Suggested change

	if ("y2" in extent) proportion /= extent.y2 - extent.y1;
	if ("x2" in extent) proportion /= extent.x2 - extent.x1;
	if ("y2" in extent && !("x2" in extent)) proportion /= extent.y2 - extent.y1;
	else if ("x2" in extent && !("y2" in extent)) proportion /= extent.x2 - extent.x1;

[Fil]

Also need to check whether value is null (like in reduceProportion), and if that is not the case, use a weighted formula.

[wirhabenzeit]

I think when binning on both x, y it can also make sense to normalise to have total area 1, see example notebook where the fill is the density which is compared to a Gaussian distribution. Without ensuring this normalisation, the density goes to zero with the bin size. For this reason I think it makes more sense to do

if ("y2" in extent) proportion /= extent.y2 - extent.y1;
if ("x2" in extent) proportion /= extent.x2 - extent.x1;

[Fil]

can you share the example?

If we want this to be consistent, then we should also normalize by hexagon area (binWidth**2) when doing hexbin?

[wirhabenzeit]

Sorry, it should be public now.

Regarding hexbin: Probably that would be consistent, although the right thing to divide by would be the area of the bin in the x-y-variable scale. However, binWidth is in pixel scale, right?

[wirhabenzeit]

Doing it consistently for hexbin would require to somehow expose the x,y-radii of the hexagonal cells. Basically, this would amount to changing hbin in src/transforms/hexbin.js to (addition highlighted with //NEW)

function hbin(data, I, X, Y, dx, sx, sy) {
    // ...
    if (bin === undefined) {
      const x = (pi + (pj & 1) / 2) * dx + ox,
        y = pj * dy + oy;
      bin = {
        index: [],
        extent: {
          data,
          x,
          y,
          rx: sx.invert(x + dx / 2) - sx.invert(x), //NEW
          ry: sy.invert(y - dy / 2) - sy.invert(y) //NEW
        }
      };
      bins.set(key, bin);
    }
    bin.index.push(i);
  }
  return bins.values();
}

and calling it with hbin(data, I, X, Y, binWidth, scales.x, scales.y).

Then the correct normalisation in reduceDensity would be

if ("rx" in extent) proportion /= 3 * extent.rx * extent.ry;

according to the area of a (hexagon)[https://en.wikipedia.org/wiki/Hexagon]

However, this adds some computational overhead by computing this for each hexagonal cell even when the values rx,ry are not needed. For linear scales rx, ry is the same for every cell but I am not sure whether the scale functions somehow expose the fact that they are linear if this is the case

[Fil]

Thanks for researching this. Note that in fact we can't guarantee anything about the x and y scales (they might not be invertible, and they might even not exist, for example when we use a projection).

This makes me think that the normalization by x2-x1 that we do here does in fact not guarantee that the scaling is correct—it will only work if x is linear (or non-existent, because identity is linear). That's the main use case we wanted to solve, but it doesn't generalize: if your histogram's base is against a log scale, the results will be inconsistent.

[wirhabenzeit]

Hmm I think the issue is that this density reduction makes sense before the scales are applied, but probably not afterwards. As far as I can see this is different between bin and hexbin.

What I mean is that for bin the density reducer as a concept also makes sense e.g. for logarithmic scales, see an example here https://observablehq.com/d/0ef7cf5eae234601 (sorry for the dots, I don't know how to produce bars in log scale). Here the reducer is applied before the scaling. The scale just influences the final appearance of the plot. Hence the reducer does not guarantee that the area under the curve is 1 after the scaling, but before the scaling. But I would argue that this is the desired behaviour because the goal is to compare the binned density with some real normalised density, irrespective of the plot scale.

[Fil]

I would argue the opposite :) If my data warrants using a log scale, it's because the thing that I am interested in analyzing is the log of the value. Scaling is not meant to make things pretty, but to connect to a deeper “truth” about the data.

This is what we do with the linear regression mark: if you use log scales, you get a linear model of the log of the values, which is equivalent to a log regression (or log-log regression if both scales are logarithmic).

To a certain extent, this is also what we do with hexbins: these hexagons are on the scaled values, and do not represent anything meaningful in data space (unless you use similar scales, for example with an aspectRatio of 1).

Linear binning is very different, as it operates in data space: the x/y bins are "rectangles", not squares. And if you want homogenous rectangles you need to manually adapt the bins so that they are equally-spaced (e.g. 1,2,4,8,16,32,64…), this is not provided by the defaults.

And, once you do that, the division by the extent of each bin should end up with the correct result? Here's what it would look like:

Plot.plot({
  marks: [
    Plot.rectY(
      pts,
      Plot.binX(
        { y: "density" },
        { x: "value", fill: "lambda", thresholds: d3.ticks(-3, 10, 20).map(d => 2**d), opacity: 0.5 }
      )
    ),
    Plot.line(density, { x: "x", y: "rho", stroke: "lambda" })
  ],
  x: { type: "log" },
  title: "Log Scale",
  grid: true,
  color: { legend: true, type: "categorical" }
})

[wirhabenzeit]

I agree that the log scale is used to reveals some truth about the data like exponential scaling.

What I meant to argue is that the reducer I am suggesting, i.e. "density of values normalised to have total area 1" should be normalised such that the area measured before the scaling and not after the scaling is equal to 1.

In my notebook I now have the same chart in three coordinate systems. Here all three charts use the "density" reducer as implemented in by pull request, i.e. with proportion /= extent.x2 - extent.x1 which is before applying the log-scale. Note that the areas look differently big, in the first plot yellow has a bigger area than blue, while in the second plot it is the other way round. However, in the space of real coordinates all areas have exactly the same size of 1, and this is what I would like the ensure.

In the applications I have in mind I plot some histogram of values (in various scalings), and compare it to the PDF of some probability measure (like exponential, Gaussian, Weibull). For this to make sense it is important that the "area under the curve" measured in the "real" coordinate system (not e.g. the log-coordinate system) is equal to 1. This is for instance also how it is done in matplotlib with plt.hist(Y, bins=100, density=True, log=True).

I understand that this becomes a bit strange for hexbin because hexbin does the binning with the scaled values other than bin which (by default) does the binning on the raw values (allowing for setting custom thresholds which can be chosen to match the scale). In order for a density reducer to make sense for hexbin there would have to be a canonical way of determining the extent of each bin in the real coordinate space.

[wirhabenzeit]

@Fil thanks for the review, will look into it! One question regarding

I think we can extend the concept to weighted density when the original channel is a value?

What do you mean by this, basically like proportion-facet but per series? What is the expected output of something like

Plot.binX({y: "density"}, {x: "xVar", y: "yVar", "stroke": "strokeVar"})

in this case? Is it the sum of the values in the bin for a given series, divided by the sum over the whole series, i.e.

sum( d.y | d.x in bin, d.z = z ) / sum( d.y | d.z = z )

or something else? I am not sure normalising by bin area makes sense in this case?

[Fil]

Normalizing with weights is for example, when a data point is a city, and represents anywhere from 1,000 to 1 million inhabitants. You want the same chart as you would have if you had one point per inhabitant.