Merge branch 'panlanfeng-master'

Author: simonbyrne

REQUIRE

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

src/KernelDensity.jl

@@ -2,12 +2,13 @@ module KernelDensity
 
 using StatsBase
 using Distributions
+using Optim
 
 import Base: conv
 import StatsBase: RealVector, RealMatrix
 import Distributions: twoπ
 
-export kde, UnivariateKDE, BivariateKDE
+export kde, UnivariateKDE, BivariateKDE, bandwidth_lscv
 
 include("univariate.jl")
 include("bivariate.jl")

src/univariate.jl

@@ -144,3 +144,36 @@ 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)-1)
+    #M is even, length(Yl2) = M/2+1 and Yl2 =[y[l]^2 for l=0 :1: M/2]
+    for j = 1:length(Yl2)-1
+        ksl = real(cf(dist,j*c))
+        zeta_star[j] = Yl2[j+1] * (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