add lscv bandwidth selector

Author: panlanfeng

REQUIRE

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

src/KernelDensity.jl

@@ -2,6 +2,7 @@ module KernelDensity
 
 using StatsBase
 using Distributions
+using Optim
 
 import Base: conv
 import StatsBase: RealVector, RealMatrix

src/univariate.jl

@@ -86,7 +86,7 @@ function tabulate(data::RealVector, midpoints::Range)
         end
     end
 
-    # returns an un-convolved KDE    
+    # returns an un-convolved KDE
     UnivariateKDE(midpoints, grid)
 end
 
@@ -126,7 +126,7 @@ end
 
 function kde(data::RealVector, dist::UnivariateDistribution;
              boundary::(Real,Real)=kde_boundary(data,std(dist)), npoints::Int=2048)
-    
+
     midpoints = kde_range(boundary,npoints)
     kde(data,midpoints,dist)
 end
@@ -144,3 +144,35 @@ function kde(data::RealVector; bandwidth=default_bandwidth(data), kernel=Normal,
     dist = kernel_dist(kernel,bandwidth)
     kde(data,dist;boundary=boundary,npoints=npoints)
 end
+
+#change the M to some larger value to get better precision of lscv
+function bandwidth_lscv(data::RealVector; kernel::DataType=Normal, M=1024)
+    n=length(data)
+    h0=default_bandwidth(data)
+    hlb = h0/sqrt(n)
+    hub = sqrt(n)*h0
+    xlb, xub = extrema(data)
+    midpoints = kde_range((xlb-4*h0, xub+4*h0), M)
+
+    k = tabulate(data, midpoints)
+    # the ft here is M/ba*sqrt(2pi) * u(s), it is M times the Yl in Silverman's book
+    Yl2 = abs2(rfft(k.density)/M)
+
+    ba = step(k.x)*M # the range b -a
+    c = -twoπ/ba
+
+    return Optim.optimize(h -> lscv(h, Yl2, kernel, c, ba, n,M), hlb, hub).minimum
+end
+
+#Silverman's book use the special case of gaussian kernel. Here the method is generalized to any symmetric kernel
+function lscv(bandwidth::Real, Yl2::RealVector, kernel::DataType, c::Real, ba::Real, n::Int,M::Int)
+    dist = kernel_dist(kernel,bandwidth)
+    zeta_star = zeros(length(Yl2))
+    for j = 1:length(Yl2)
+        ksl = real(cf(dist,j*c))
+        zeta_star[j] = Yl2[j] * (ksl * ksl - 2 * ksl)
+    end
+    #Correct the error in silverman's book
+    #∫ (cf^2 -2cf)u(s)²ds <- ∑(cf^2 - 2cf)*Yl2*ba²/2pi * c
+    sum(zeta_star) * abs(c)*ba*ba/(2*pi) + pdf(dist, 0.0)/n
+end