Density estimates: geom_density

Mix.install([
  {:giggly, path: Path.join(__DIR__, "..")},
  {:kino, "~> 0.19"}
  ])

alias Giggly.Plot

diamonds = Giggly.Datasets.diamonds()

render = fn plot ->
  plot
  |> Giggly.to_svg()
  |> IO.iodata_to_binary()
  |> Kino.Image.new(:svg)
end

About this notebook

This notebook follows the examples on ggplot2's reference page for geom_density(), which covers stat_density() too, in the order the page gives them. Each example shows ggplot2's code first and then Giggly's. An example Giggly cannot draw yet keeps its ggplot2 code and names what it is waiting for, so the page is also a list of what remains.

geom_density() draws a kernel density estimate, a smoothed version of the histogram, which suits continuous data drawn from a smooth distribution; stat_density() computes it. Besides the density, the stat computes counts, the density times the number of rows, for stacking densities. Giggly has both, Giggly.Plot.geom_density/3 and Giggly.Plot.stat_density/3, which take ggplot2's options as keyword options (bw:, adjust:, kernel:, n:, trim:, bounds: and orientation:), stack curves with position: :stack, and map what the stat computes, such as the count, with {:after_stat, [:count]}. It has no position_fill(), and of the bandwidth rules bw names, only R's "nrd0" and "nrd", as :nrd0 and :nrd, not "ucv", "bcv" or "SJ". A plot is described with Giggly.Plot and drawn with Giggly.to_svg/2, which returns an SVG document. The setup cell aliases Giggly.Plot, so the cells write Plot.new/2, and defines render, which draws a plot and shows it below the cell with Kino.Image. Every example uses ggplot2's diamonds, Giggly.Datasets.diamonds/0, which the setup cell binds to diamonds, ggplot2's name for it.

A density curve

ggplot(diamonds, aes(carat)) +
  geom_density()

diamonds
|> Plot.new(%{x: :carat})
|> Plot.geom_density()
|> render.()

Horizontal densities

# Map the values to y to flip the orientation
ggplot(diamonds, aes(y = carat)) +
  geom_density()

diamonds
|> Plot.new(%{y: :carat})
|> Plot.geom_density()
|> render.()

Rougher and smoother curves

ggplot(diamonds, aes(carat)) +
  geom_density(adjust = 1/5)
ggplot(diamonds, aes(carat)) +
  geom_density(adjust = 5)

adjust multiplies the bandwidth the density is smoothed over.

diamonds
|> Plot.new(%{x: :carat})
|> Plot.geom_density(adjust: 1 / 5)
|> render.()

diamonds
|> Plot.new(%{x: :carat})
|> Plot.geom_density(adjust: 5)
|> render.()

A curve for each cut

ggplot(diamonds, aes(depth, colour = cut)) +
  geom_density() +
  xlim(55, 70)
ggplot(diamonds, aes(depth, fill = cut, colour = cut)) +
  geom_density(alpha = 0.1) +
  xlim(55, 70)

The limits drop the diamonds whose depth lies outside them before each curve is estimated.

diamonds
|> Plot.new(%{x: :depth, colour: :cut})
|> Plot.geom_density()
|> Plot.xlim(55, 70)
|> render.()

diamonds
|> Plot.new(%{x: :depth, fill: :cut, colour: :cut})
|> Plot.geom_density(alpha: 0.1)
|> Plot.xlim(55, 70)
|> render.()

Bounded data

# Use `bounds` to adjust computation for known data limits
big_diamonds <- diamonds[diamonds$carat >= 1, ]
ggplot(big_diamonds, aes(carat)) +
  geom_density(color = 'red') +
  geom_density(bounds = c(1, Inf), color = 'blue')

The diamonds of a carat or more have no values below 1, which the red curve, knowing no bounds, smooths past. The cell keeps the rows R's indexing keeps, each column's values where carat is 1 or more, and writes R's Inf as :infinity.

kept = Enum.map(Keyword.fetch!(diamonds, :carat), fn carat -> carat >= 1 end)

big_diamonds =
  Enum.map(diamonds, fn {name, values} ->
    rows =
      values
      |> Enum.zip(kept)
      |> Enum.filter(fn {_value, keep} -> keep end)
      |> Enum.map(fn {value, _keep} -> value end)

    {name, rows}
  end)

big_diamonds
|> Plot.new(%{x: :carat})
|> Plot.geom_density(colour: "red")
|> Plot.geom_density(bounds: {1, :infinity}, colour: "blue")
|> render.()

Stacked densities

# Stacked density plots: if you want to create a stacked density plot, you
# probably want to 'count' (density * n) variable instead of the default
# density
# Loses marginal densities
ggplot(diamonds, aes(carat, fill = cut)) +
  geom_density(position = "stack")
# Preserves marginal densities
ggplot(diamonds, aes(carat, after_stat(count), fill = cut)) +
  geom_density(position = "stack")

diamonds
|> Plot.new(%{x: :carat, fill: :cut})
|> Plot.geom_density(position: :stack)
|> render.()

diamonds
|> Plot.new(%{x: :carat, y: {:after_stat, [:count]}, fill: :cut})
|> Plot.geom_density(position: :stack)
|> render.()

Conditional densities

# You can use position="fill" to produce a conditional density estimate
ggplot(diamonds, aes(carat, after_stat(count), fill = cut)) +
  geom_density(position = "fill")

Not in Giggly yet: position_fill().