# R function bw.ss for a new bandwidth selector for kernel density #An example in the bottom of this file ######################################################## U.func = function(h, s, q, w, mu, sigma2) { k = length(w) sigma.local = outer(sigma2, sigma2, "+") + q*h*h sigma.local = sqrt(sigma.local) diff1 = outer(mu, mu, "-")/sigma.local term1 = K.deriv(diff1, 2*s)/sigma.local^(2*s+1) (1/k^2)*(h^(2*s))*sum(term1) } MISE.cal.exact = function(h, n, w, mu, sigma2) { term1 = 1/(2*sqrt(pi))/(n*h) term2 = (1-1/n)*U.func(h, 0, 2, w, mu, sigma2) term3 = -2*U.func(h, 0, 1, w, mu, sigma2) term4 = U.func(h, 0, 0, w, mu, sigma2) term1 + term2 + term3 + term4 } V.deriv.func = function(h, s, q, w, mu, sigma2) { k = length(w) sigma.local = outer(sigma2, sigma2, "+") + q*h*h sigma.local = sqrt(sigma.local) diff1 = outer(mu, mu, "-")/sigma.local term1 = K.deriv(diff1, 2*s)/(sigma.local^(2*s+1)) term2 = K.deriv(diff1, 2*s+2)/(sigma.local^(2*s+3)) term1 = 2*s*h^(2*s-1)*sum(term1) term2 = q*h^(2*s+1)*sum(term2) (1/k^2)*(term1+term2) } MISE.deriv.cal.exact = function(h, n, w, mu, sigma2) { term1 = -1/(2*sqrt(pi))/(n*h*h) term2 = (1-1/n)*V.deriv.func(h, 0, 2, w, mu, sigma2) term3 = -2*V.deriv.func(h, 0, 1, w, mu, sigma2) term1 + term2 + term3 } #################################### H = function(u, k) { if(k==0) 1 else if(k==1) u else if(k==2) u^2 - 1 else if(k==4) u^4 - 6*u^2 + 3 else if(k==6) u^6 - 15*u^4 + 45*u^2 - 15 else if(k==8) u^8 - 28*u^6 + 210*u^4 - 420*u^2 + 105 else if(k==10) u^10 - 45*u^8 + 630*u^6 - 3150*u^4 + 4725*u^2 - 945 else stop("invalid k") } K.deriv = function(u, k) H(u, k)*dnorm(u) Phi = function(x, h1, k) { n = length(x) x = x/h1 temp = outer(x, x, "-") n/(n-1)*mean(K.deriv(temp, k))/h1^(k+1) } Phi.population.normal = function(k, sigma) (-1)^(k/2)*factorial(k)/(2*sigma)^(k+1)/factorial(k/2)/sqrt(pi) g.optimal = function(Phi.population.2higer, k, n) (-2*K.deriv(0, k)/Phi.population.2higer/n)^(1/(k+3)) ##################################################################### bw.liao.exact = function(x) { n = length(x) RK = 1/2/sqrt(pi) temp = 2*K.deriv(0, 4)/RK g4.initial = function(x) { h1 = bw.SJ(x) func1 = function(log.g4) { g4 = exp(log.g4) g6 = .992*g4*n^(1/7-1/9) #for the ratio of Phi4 over Phi6 ratio = -Phi(x, g6, 4)/Phi(x, g6, 6) theFactor = (temp*ratio)^(1/7) (g4/theFactor)^(7/5) - h1 } interval = log(c(h1/10, h1*10)) exp(uniroot(func1, interval)$root) } h1 = NULL func1 = function(log.g4) #g4 is the bandwidth to estimate R2 from sample x_1,...,x_n { g4 = exp(log.g4) g6 = .992*g4*n^(1/7-1/9) ratio = -Phi(x, g6, 4)/Phi(x, g6, 6) theFactor = (temp*ratio)^(1/7) #reference distribution with R2 h1 <<- (g4/theFactor)^(7/5) g0 = g4/sqrt(2) sigma2 = rep(g0^2, n) #for using Wand and Marron's exact formula in Theorem 5.1 w = rep(1/n, n) MISE.deriv.cal.exact(h1, n, w, x, sigma2) } g4 = g4.initial(x) interval = log(c(g4/5, g4*5)) value1 = func1(interval[1]) value2 = func1(interval[2]) if(value1*value2 > 0) bw.SJ(x) else { uniroot(func1, interval, f.lower=value1, f.upper=value2, tol = 10*.Machine$double.eps^0.25) h1 } } #### the above is exact solution and below is approximate solution ######### bw.liao.approx = function (x, nb = 4000, SJ=F) #nb is the number of bins { n = length(x) SDh <- function(x, h, n, d) .C(stats:::R_band_phi4_bin, as.integer(n), #Phi4 x will be replaced by a processed version cnt as.integer(length(x)), as.double(d), x, as.double(h), u = double(1))$u TDh <- function(x, h, n, d) .C(stats:::R_band_phi6_bin, as.integer(n), #Phi6 x will be replaced by a processed version cnt as.integer(length(x)), as.double(d), x, as.double(h), u = double(1))$u Z <- .C(stats:::R_band_den_bin, as.integer(n), as.integer(nb), d = double(1), x, cnt = integer(nb)) d <- Z$d #needed for the faster functions cnt <- as.integer(Z$cnt) #this will be a replacement for x in SDh and TDh RK = 1/2/sqrt(pi) temp = 2*K.deriv(0, 4)/RK #### Below is for SJ selector ################ #hmax <- 1.144 * sqrt(var(x)) * n^(-1/5) #scale <- min(sqrt(var(x)), IQR(x)/1.349) #g4 <- 1.24 * scale * n^(-1/7) #g6 <- 1.23 * scale * n^(-1/9) #ratio0 = -SDh(cnt, g4, n, d)/TDh(cnt, g6, n, d) ##### Above is for SJ selector############################################## ########################################## g4.initial = function(x) { h0 = bw.SJ(x) func1 = function(log.g4) { g4 = exp(log.g4) g6 = .992*g4*n^(1/7-1/9) #for the ratio of Phi4 over Phi6 ratio = -SDh(cnt, g6, n, d)/TDh(cnt, g6, n, d) theFactor = (temp*ratio)^(1/7) (g4/theFactor)^(7/5) - h0 } interval = log(c(h0/10, h0*10)) exp(uniroot(func1, interval)$root) } h1 = NULL key.func1 = function(log.g4) { g4 = exp(log.g4) g6 = .992*g4*n^(1/7-1/9) #theFactor0 = (temp*ratio0)^(1/7) #SJ #h0.0 = (g4/theFactor0)^(7/5) #SJ ratio = -SDh(cnt, g6, n, d)/TDh(cnt, g6, n, d) theFactor = (temp*ratio)^(1/7) h0.1 = (g4/theFactor)^(7/5) h0.2 = (RK/SDh(cnt, g4, n, d)/n)^.2 ratio2 = -SDh(cnt, g4, n, d)/TDh(cnt, g4, n, d) J2 = 3/20/ratio2 h0.3 = h0.2 + J2*h0.2^3 # h1 <<- h0.0; h0.0 - h0.2 #SJ selector h1 <<- h0.1 h0.1 - h0.3 #Liao selector } g4 = g4.initial(x) interval = log(c(g4/5, g4*5)) value1 = key.func1(interval[1]) value2 = key.func1(interval[2]) if(value1*value2 > 0) bw.SJ(x) else { uniroot(key.func1, interval, f.lower=value1, f.upper=value2) h1 } } ###########################################################################3 # Liao, Wu and Lin (2009) Improving Sheather and Jones bandwidth selector in kernel density estimation bw.liao = function(x) { h0 = bw.SJ(x) if(length(x) <= 400) h2.value = bw.liao.exact(x) else h2.value = bw.liao.approx(x) h3.value = min(h2.value, h0) print("SJ bandwidth versus our bandwidth") print(c(h0, h3.value)) h3.value } ################################################# ###### Example #################### x = rnorm(200) bw.liao(x)