Add interpolation functionality, closes #9

Author: simonbyrne

README.md

@@ -9,21 +9,35 @@ Kernel density estimators for julia.
 
 ### Univariate
 The main accessor function is `kde`:
+
 ```
 kde(data)
 ```
-will construct a `UnivariateKDE` object from the real vector `data`. The optional keyword arguments are
+
+will construct a `UnivariateKDE` object from the real vector `data`. The
+optional keyword arguments are
 * `boundary`: the lower and upper limits of the kde as a tuple. Due to the
   fourier transforms used internally, there should be sufficient spacing to
   prevent wrap-around at the boundaries.
 * `npoints`: the number of interpolation points to use. The function uses
   fast Fourier transforms (FFTs) internally, so for optimal efficiency this
   should be a power of 2 (default = 2048).
-* `kernel`: the distributional family to use as the kernel (default =
-  `Normal`). To add your own kernel, extend the internal `kernel_dist` function.
+* `kernel`: the distributional family from
+  [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) to use as
+  the kernel (default = `Normal`). To add your own kernel, extend the internal
+  `kernel_dist` function.
 * `bandwidth`: the bandwidth of the kernel. Default is to use Silverman's
   rule.
 
+A related function
+
+``` kde_lscv(data) ```
+
+will construct a `UnivariateKDE` object, with the bandwidth selected by
+least-squares cross validation. It accepts the above keyword arguments, except
+`bandwidth`.
+
+
 There are also some slightly more advanced interfaces:
 ```
 kde(data, midpoints::Range)
@@ -56,6 +70,35 @@ kde(datamatrix)
 Similarly, the optional arguments all now take tuple arguments:
 e.g. `boundary` now takes a tuple of tuples `((xlo,xhi),(ylo,yhi))`.
 
+### Interpolation
+
+The KDE objects are stored as gridded density values, with attached
+coordinates. These are typically sufficient for plotting (see below), but
+intermediate values can be interpolated using the
+[Grid.jl](https://github.com/timholy/Grid.jl) package via the `pdf` method
+(extended from Distributions.jl).
+
+```
+pdf(k::UnivariateKDE, x)
+pdf(k::BivariateKDE, x, y)
+```
+
+where `x` and `y` are real numbers or arrays.
+
+If you are making multiple calls to `pdf`, it will be more efficient to
+construct an intermediate `InterpKDE` to store the intermediate `InterpGrid`
+object:
+
+```
+ik = InterpKDE(k)
+pdf(ik, x)
+```
+
+`InterpKDE` can also take additional arguments specifying the
+`BoundaryCondition` (default=`BCnan`) and `InterpType` (default=
+`InterpQuadratic`).
+
+
 ## Plotting
 
 The [Winston.jl](https://github.com/nolta/Winston.jl) and

REQUIRE

@@ -2,3 +2,4 @@ julia 0.3-
 StatsBase
 Distributions 0.4.6
 Optim
+Grid

src/KernelDensity.jl

@@ -3,15 +3,17 @@ module KernelDensity
 using StatsBase
 using Distributions
 using Optim
+using Grid
 
 import Base: conv
 import StatsBase: RealVector, RealMatrix
-import Distributions: twoπ
+import Distributions: twoπ, pdf
 
-export kde, kde_lscv, UnivariateKDE, BivariateKDE
+export kde, kde_lscv, UnivariateKDE, BivariateKDE, InterpKDE, pdf
 
 include("univariate.jl")
 include("bivariate.jl")
+include("interp.jl")
 
 macro glue(pkg)
     path = joinpath(dirname(Base.source_path(nothing)),"glue",string(pkg,".jl"))

src/bivariate.jl

@@ -1,5 +1,5 @@
 # Store both grid and density for KDE over R2
-immutable BivariateKDE{Rx<:Range,Ry<:Range}
+type BivariateKDE{Rx<:Range,Ry<:Range}
     x::Rx
     y::Ry
     density::Matrix{Float64}

src/interp.jl

@@ -0,0 +1,19 @@
+type InterpKDE{K,G}
+    kde::K
+    grid::G
+end
+
+function InterpKDE{BC<:BoundaryCondition, IT<:InterpType}(k::UnivariateKDE, bc::Type{BC}=BCnan,it::Type{IT}=InterpQuadratic)
+    g = CoordInterpGrid(k.x, k.density, bc, it)
+    InterpKDE(k,g)
+end
+function InterpKDE{BC<:BoundaryCondition, IT<:InterpType}(k::BivariateKDE, bc::Type{BC}=BCnan,it::Type{IT}=InterpQuadratic)
+    g = CoordInterpGrid((k.x,k.y),k.density,bc,it)
+    InterpKDE{typeof(k),typeof(g)}(k,g)
+end
+
+
+pdf(ik::InterpKDE,x...) = ik.grid[x...]
+
+pdf(k::UnivariateKDE,x) = pdf(InterpKDE(k),x)
+pdf(k::BivariateKDE,x,y) = pdf(InterpKDE(k),x,y)

src/univariate.jl

@@ -1,5 +1,5 @@
 # Store both grid and density for KDE over the real line
-immutable UnivariateKDE{R<:Range}
+type UnivariateKDE{R<:Range}
     x::R
     density::Vector{Float64}
 end

test/interp.jl

@@ -0,0 +1,12 @@
+using Base.Test
+using KernelDensity
+using Grid
+
+X = randn(100)
+Y = randn(100)
+
+k = kde(X)
+@test_approx_eq pdf(k, k.x) k.density
+
+k = kde((X,Y))
+@test_approx_eq pdf(k, k.x, k.y) k.density