Merge pull request #1103 from Mattriks/Geom_ellipse

Geom.ellipse

.travis.yml

@@ -2,7 +2,6 @@ language: julia
 os:
   - linux
 julia:
-  - 0.5
   - 0.6
   - nightly
 notifications:

docs/src/lib/geoms/geom_ellipse.md

@@ -0,0 +1,47 @@
+```@meta
+Author = "Mattriks"
+```
+
+# Geom.ellipse
+
+Confidence ellipse for a scatter or group of points, using a parametric multivariate distribution e.g. multivariate normal. `Geom.ellipse` is an instance of [`Geom.polygon`](@ref)
+
+## Aesthetics
+
+  * `x`: Position of points.
+  * `y`: Position of points.
+  * `color` (optional): Color.
+  * `group` (optional): Group.
+
+## Arguments
+
+  * `distribution`: A multivariate distribution. Default is `MvNormal`.
+  * `levels`: The quantiles for which confidence ellipses are calculated. Default is [0.95].
+  * `nsegments`: Number of segments to draw each ellipse. Default is 51.
+
+
+## Examples
+
+```@setup 1
+using RDatasets, Gadfly
+Gadfly.set_default_plot_size(14cm, 8cm)
+```
+
+```@example 1
+D = dataset("datasets","faithful")
+D[:g] = D[:Eruptions].>3.0
+
+coord = Coord.cartesian(ymin=35, ymax=100)
+
+pa = plot(D, coord,
+    x=:Eruptions, y=:Waiting, group=:g,
+    Geom.point, Geom.ellipse
+)
+pb = plot(D, coord,
+    x=:Eruptions, y=:Waiting, color=:g,
+    Geom.point, Geom.ellipse,
+    layer(Geom.ellipse(levels=[0.99]), style(line_style=:dot)),
+    style(key_position=:none), Guide.ylabel(nothing)
+)
+hstack(pa,pb)
+```

src/geom/polygon.jl

@@ -1,17 +1,25 @@
 immutable PolygonGeometry <: Gadfly.GeometryElement
+    default_statistic::Gadfly.StatisticElement
     order::Int
     fill::Bool
     preserve_order::Bool
     tag::Symbol
 end
 
-PolygonGeometry(; order=0, fill=false, preserve_order=false, tag=empty_tag) =
-        PolygonGeometry(order, fill, preserve_order, tag)
+PolygonGeometry(default_statistic=Gadfly.Stat.identity(); order=0, fill=false, preserve_order=false, tag=empty_tag) =
+        PolygonGeometry(default_statistic, order, fill, preserve_order, tag)
 
 const polygon = PolygonGeometry
 
 element_aesthetics(::PolygonGeometry) = [:x, :y, :color, :group]
 
+ellipse(;distribution::(@compat Type{<:ContinuousMultivariateDistribution})=MvNormal,
+    levels::Vector=[0.95], nsegments::Int=51, fill::Bool=false) =
+    PolygonGeometry(Gadfly.Stat.ellipse(distribution, levels, nsegments), preserve_order=true, fill=fill)
+
+default_statistic(geom::PolygonGeometry) = geom.default_statistic
+
+
 function polygon_points(xs, ys, preserve_order)
     T = (Tuple{eltype(xs), eltype(ys)})
     if preserve_order
@@ -36,6 +44,8 @@ function render(geom::PolygonGeometry, theme::Gadfly.Theme,
     ctx = context(order=geom.order)
     T = (eltype(aes.x), eltype(aes.y))
 
+    line_style = Gadfly.get_stroke_vector(theme.line_style)
+
     if aes.group != nothing
         XT, YT = eltype(aes.x), eltype(aes.y)
         xs = DefaultDict{Any, Vector{XT}}(() -> XT[])
@@ -86,5 +96,5 @@ function render(geom::PolygonGeometry, theme::Gadfly.Theme,
         end
     end
 
-    return compose!(ctx, linewidth(theme.line_width), svgclass("geometry"))
+    return compose!(ctx, linewidth(theme.line_width), strokedash(line_style), svgclass("geometry"))
 end

src/geometry.jl

@@ -5,6 +5,7 @@ using Compat
 using Compose
 using DataArrays
 using DataStructures
+using Distributions
 using Gadfly
 using Measures
 

src/statistics.jl

@@ -8,14 +8,15 @@ using Compat
 using Compose
 using DataArrays
 using DataStructures
+using Distributions
 using Hexagons
 using Loess
 using CoupledFields # It is registered in METADATA.jl
 
 import Gadfly: Scale, Coord, input_aesthetics, output_aesthetics,
                default_scales, isconcrete, nonzero_length, setfield!
 import KernelDensity
-import Distributions: Uniform, Distribution, qqbuild
+# import Distributions: Uniform, Distribution, qqbuild
 import IterTools: chain, distinct
 import Compat.Iterators: cycle, product
 
@@ -1712,6 +1713,7 @@ function apply_statistic(stat::EnumerateStatistic,
     end
 end
 
+### Vector Field Statistic 
 
 immutable VecFieldStatistic <: Gadfly.StatisticElement
     smoothness::Float64
@@ -1779,6 +1781,7 @@ function apply_statistic(stat::VecFieldStatistic,
 
 end
 
+### Hair Statistic
 
 immutable HairStatistic <: Gadfly.StatisticElement
     intercept
@@ -1808,6 +1811,74 @@ function apply_statistic(stat::HairStatistic,
 end
 
 
+### Ellipse Statistic
+
+immutable EllipseStatistic <: Gadfly.StatisticElement
+    distribution::@compat Type{<:ContinuousMultivariateDistribution}
+    levels::@compat Vector{<:AbstractFloat}
+    nsegments::Int
+end
+
+function EllipseStatistic(; 
+        distribution::(@compat Type{<:ContinuousMultivariateDistribution})=MvNormal,
+        levels::Vector{Float64}=[0.95], 
+        nsegments::Int=51 )
+    return EllipseStatistic(distribution, levels, nsegments)
+end
+
+Gadfly.input_aesthetics(stat::EllipseStatistic) = [:x, :y]
+Gadfly.output_aesthetics(stat::EllipseStatistic) = [:x, :y]
+Gadfly.default_scales(stat::EllipseStatistic) = [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous()]
+
+const ellipse = EllipseStatistic
+
+function Gadfly.Stat.apply_statistic(stat::EllipseStatistic,
+                         scales::Dict{Symbol, Gadfly.ScaleElement},
+                         coord::Gadfly.CoordinateElement,
+                         aes::Gadfly.Aesthetics)
+
+    Dat = [aes.x aes.y]
+    grouped_xy = Dict(1=>Dat)
+    grouped_color = Dict{Int, Gadfly.ColorOrNothing}(1=>nothing)
+    colorflag = aes.color != nothing
+    aes.group = (colorflag ? aes.color : aes.group)
+
+    if aes.group != nothing
+        ug = unique(aes.group) 
+        grouped_xy = Dict(g=>Dat[aes.group.==g,:] for g in ug)
+        grouped_color = Dict(g=>first(aes.group[aes.group.==g]) for g in ug)
+    end
+
+    levels = Float64[]
+    colors = eltype(aes.color)[]
+    ellipse_x = eltype(Dat)[]
+    ellipse_y = eltype(Dat)[]
+
+    dfn = 2
+    θ = 2π*(0:stat.nsegments)/stat.nsegments
+    n = length(θ)
+    for (g, data) in grouped_xy
+        dfd = size(data,1)-1
+        dhat = fit(stat.distribution, data')
+        Σ½ = chol(cov(dhat))
+        rv = sqrt.(dfn*[quantile(FDist(dfn,dfd), p) for p in stat.levels])
+        ellxy =  [cos.(θ) sin.(θ)] * Σ½
+        μ = mean(dhat)
+        for r in rv
+            append!(ellipse_x, r*ellxy[:,1].+μ[1])
+            append!(ellipse_y, r*ellxy[:,2].+μ[2])
+            append!(colors, fill(grouped_color[g], n))
+            append!(levels, fill(r, n))
+        end
+    end
+
+    aes.group = PooledDataArray(levels)
+    colorflag && (aes.color = colors)
+    aes.x = ellipse_x
+    aes.y = ellipse_y
+end
+
+
 
 
 end # module Stat

test/testscripts/ellipse.jl

@@ -0,0 +1,9 @@
+using Distributions, Gadfly
+set_default_plot_size(6.6inch, 3.3inch)
+
+srand(123)
+d = rand(MvNormal([2, 2],[1.0 0.7; 0.7 1.0]), 50)'
+
+pa= plot(x=d[:,1], y=d[:,2], Geom.point, layer(Stat.ellipse, Geom.polygon(preserve_order=true)))
+pb= plot(x=d[:,1], y=d[:,2], Geom.point, Geom.ellipse(levels=[0.95, 0.99]))
+hstack(pa,pb)