Quantile bars

NEWS.md

@@ -3,7 +3,7 @@ Each release typically has a number of minor bug fixes beyond what is listed her
 
 # Version 1.x
 
-
+ * Add `Stat.quantile_bars` (#1521)
  * Add "vertical" orientation for `Geom.ribbon` (#1513)
  * Support one-length aesthetics for `Geom.polygon` and `Geom.ribbon` (#1511)
  * Enable `color` grouping for `Geom.density2d` (#1508)

docs/Project.toml

@@ -1,14 +1,15 @@
 [deps]
 Cairo = "159f3aea-2a34-519c-b102-8c37f9878175"
-Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004"
 RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
 Showoff = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
 
 [compat]
+CategoricalArrays = "0.9"
 Documenter = "0.26"

docs/make.jl

@@ -3,7 +3,8 @@ using Documenter, Gadfly, Compose, Cairo
 makedocs(
     modules = [Gadfly],
     format = Documenter.HTML(
-        assets = ["assets/favicon.ico"]
+        assets = ["assets/favicon.ico"],
+        prettyurls = get(ENV, "CI", nothing) == "true"
     ),
     clean = false,
     sitename = "Gadfly.jl",

docs/src/gallery/statistics.md

@@ -33,6 +33,24 @@ p2 = plot(vcat(Db...), x=:x, color=:u, Theme(alphas=[0.6]),
 hstack(p1,p2)
 ```
 
+## [`Stat.quantile_bars`](@ref)
+
+```@example
+using CategoricalArrays
+using Gadfly
+set_default_plot_size(14cm, 8cm)
+n = 400
+group = repeat([-1, 1], inner=200)
+x = randn(n) .+ group
+
+plot(x=x, color=categorical(group), Guide.colorkey(title="", pos=[3.6,0.7]),
+    layer(Stat.density, Geom.line, Geom.polygon(fill=true, preserve_order=true), alpha=[0.4]),
+    layer(Stat.quantile_bars(quantiles=[0.05, 0.95]), Geom.segment),
+    Guide.title("Density with bars showing the central 90% CI"),
+    Guide.ylabel("Density"), Coord.cartesian(xmin=-4, xmax=4)
+)
+```
+
 ## [`Stat.dodge`](@ref)
 
 ```@example

src/statistics.jl

@@ -2144,9 +2144,83 @@ function apply_statistic(stat::DodgeStatistic,
     end
 end
 
+struct QuantileBarsStatistic <: Gadfly.StatisticElement
+    quantiles::Vector{Float64}
+    # We cannot avoid these by combining our statistic with Stat.density,
+    # because we need the raw data as well as the kernel density.
+    n::Int # Number of points sampled.
+    bw::Real # Bandwidth used for the kernel density estimation.
+end
+QuantileBarsStatistic(; quantiles=[0.025, 0.975], n=256, bandwidth=-Inf) =
+    QuantileBarsStatistic(quantiles, n, bandwidth)
+
+input_aesthetics(stat::QuantileBarsStatistic) = [:x]
+output_aesthetics(stat::QuantileBarsStatistic) = [:x, :y, :xend, :yend]
+
+"""
+    Stat.quantile_bars[(; quantiles=[0.025, 0.975], bar_width=0.1, n=256, bandwidth=-Inf)]
+
+Transform the point in $(aes2str(input_aesthetics(quantile_bars()))) into a set of
+$(aes2str(output_aesthetics(quantile_bars()))) points. These points can then be drawn
+via [`Geom.segment`](@ref Gadfly.Geom.segment). Here, `bandwidth` works independently
+from the `bandwidth` setting for `Stat.density`.
+"""
+const quantile_bars = QuantileBarsStatistic
+
+"""
+    _calculate_quantile_bar(stat::QuantileBarsStatistic, aes)
 
+Helper function for `apply_statistic(stat::QuantileBarsStatistic, ...)`.
+"""
+function _calculate_quantile_bar(stat::QuantileBarsStatistic, xs)
+    isa(xs[1], Real) || error("Kernel density estimation only works on Real types.")
+
+    window = stat.bw <= 0.0 ? KernelDensity.default_bandwidth(xs) : stat.bw
+    k = KernelDensity.kde(xs, bandwidth=window, npoints=stat.n)
 
+    x = quantile(xs, stat.quantiles)
+    y = zeros(length(x))
+    xend = x
+    yend = pdf(k, x)
+
+    return x, y, xend, yend
+end
 
+function apply_statistic(stat::QuantileBarsStatistic,
+                       scales::Dict{Symbol, Gadfly.ScaleElement},
+                       coord::Gadfly.CoordinateElement,
+                       aes::Gadfly.Aesthetics)
+    Gadfly.assert_aesthetics_defined("QuantileBarsStatistic", aes, :x)
 
+    if aes.color === nothing
+        aes.x, aes.y, aes.xend, aes.yend = _calculate_quantile_bar(stat, aes.x)
+    else
+        groups = Dict()
+        for (x, c) in zip(aes.x, Gadfly.cycle(aes.color))
+            if !haskey(groups, c)
+                groups[c] = Float64[x]
+            else
+                push!(groups[c], x)
+            end
+        end
+
+        colors = Array{Gadfly.RGB{Float32}}(undef, 0)
+        aes.x = Array{Float64}(undef, 0)
+        aes.y = Array{Float64}(undef, 0)
+        aes.xend = Array{Float64}(undef, 0)
+        aes.yend = Array{Float64}(undef, 0)
+        for (c, xs) in groups
+            x, y, xend, yend = _calculate_quantile_bar(stat, xs)
+
+            append!(aes.x, x)
+            append!(aes.y, y)
+            append!(aes.xend, xend)
+            append!(aes.yend, yend)
+            append!(colors, fill(c, length(x)))
+        end
+        aes.color = discretize_make_ia(colors)
+    end
+    aes.y_label = Gadfly.Scale.identity_formatter
+end
 
 end # module Stat

test/testscripts/stat_quantile_bars.jl

@@ -0,0 +1,15 @@
+using Gadfly, RDatasets, Distributions
+
+set_default_plot_size(6inch, 3inch)
+
+xs = dataset("datasets", "faithful").Waiting
+
+quantiles = [0.1, 0.9]
+p1 = plot(x=xs, Geom.density, Guide.title("quantiles=$quantiles"),
+        layer(x=xs, Stat.quantile_bars(quantiles=quantiles), Geom.segment))
+
+bandwidth = 0.6
+p2 = plot(x=xs, Stat.density(bandwidth=bandwidth), Geom.line, Guide.title("bandwidth=$bandwidth"),
+    layer(x=xs, Stat.quantile_bars(bandwidth=bandwidth), Geom.segment))
+
+p = hstack(p1, p2)