Nothing
###############################################################################
## Estimate g_AMSE pilot bandwidths for even orders - 2-dim
##
## Parameters
## r - (r1, r2) partial derivative
## n - sample size
## psi1 - psi_(r + (2,0))
## psi2 - psi_(r + (0,2))
##
## Returns
## g_AMSE pilot bandwidths for even orders
###############################################################################
gamse.even.2d <- function(r, n, psi1, psi2)
{
d <- 2
num <- -2 * dmvnorm.deriv(x=c(0,0), deriv.order=r, Sigma=diag(2), deriv.vec=FALSE)
den <- (psi1 + psi2) * n
g.amse <- (num/den)^(1/(2 + d + sum(r)))
return(g.amse)
}
###############################################################################
## Estimate g_AMSE pilot bandwidths for odd orders - 2-dim
##
## Parameters
##
## r - (r1, r2) partial derivative
## n - sample size
## psi1 - psi_(r + (2,0))
## psi2 - psi_(r + (0,2))
## psi00 - psi_(0,0)
## RK - R(K^(r))
##
## Returns
## g_AMSE pilot bandwidths for odd orders
###############################################################################
gamse.odd.2d <- function(r, n, psi1, psi2, psi00, RK)
{
d <- 2
num <- 2 * psi00 * (2 * sum(r) + d) * RK
den <- (psi1 + psi2)^2 * n^2
g.amse <- (num/den)^(1/(2*sum(r) + d + 4))
return(g.amse)
}
###############################################################################
## Estimate g_SAMSE pilot bandwidth - 2- to 6-dim
##
## Parameters
## Sigma.star - scaled variance matrix
## n - sample size
##
## Returns
## g_SAMSE pilot bandwidth
###############################################################################
gsamse <- function(Sigma.star, n, modr, nstage=1, psihat=NULL)
{
d <- ncol(Sigma.star)
K <- numeric(); psi <- numeric()
## 4th order g_SAMSE
K <- dmvnorm.deriv(x=rep(0,d), deriv.order=modr, Sigma=diag(d), add.index=TRUE, deriv.vec=FALSE)
K <- K$deriv[apply(K$deriv.ind, 1, is.even)]
if (modr==4)
{
derivt4 <- dmvnorm.deriv(x=rep(0,d), deriv.order=4, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
derivt6 <- dmvnorm.deriv(x=rep(0,d), deriv.order=6, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
for (i in 1:nrow(derivt4))
{
r <- derivt4[i,]
if (is.even(r))
{
A3psi <- 0
for (j in 1:d)
{
if (nstage==1) A3psi <- A3psi + psins(r=r+2*elem(j,d), Sigma=Sigma.star)
else if (nstage==2) A3psi <- A3psi + psihat[which.mat(r=r+2*elem(j,d), mat=derivt6)]
}
psi <- c(psi, A3psi)
}
}
}
## 6th order g_SAMSE
else if (modr==6)
{
derivt6 <- dmvnorm.deriv(x=rep(0,d), deriv.order=6, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
for (i in 1:nrow(derivt6))
{
r <- derivt6[i,]
if (is.even(r))
{
A3psi <- 0
for (j in 1:d)
A3psi <- A3psi + psins(r=r+2*elem(j,d), Sigma=Sigma.star)
psi <- c(psi, A3psi)
}
}
}
## see thesis for formula
A1 <- sum(K^2)
A2 <- sum(K * psi)
A3 <- sum(psi^2)
B1 <- (2*modr + 2*d)*A1
B2 <- (modr + d - 2)*A2
B3 <- A3
gamma1 <- (-B2 + sqrt(B2^2 + 4*B1*B3)) / (2*B1)
g.samse <- (gamma1 * n)^(-1/(modr + d + 2))
return (g.samse)
}
##############################################################################
## Scalar pilot selector for derivatives r>0 from Chacon & Duong (2011)
## Generalisation of gsamse for r>0
##############################################################################
gdscalar <- function(x, d, r, n, verbose, binned=FALSE, nstage=1, scv=FALSE)
{
if (scv) cf <- c(2^(-d), 2^(-d/2+1), 4)
else cf <- c(1,1,1)
S <- var(x)
Sinv <- chol2inv(chol(S))
if (nstage==1)
{
psi2r6.ns <- psins(r=2*r+6, Sigma=S)
B3 <- norm(Matrix(psi2r6.ns, nrow=1) %*% (Matrix(vec(diag(d)), ncol=1) %*% Matrix(vec(diag(d)), nrow=1) %x% Diagonal(d^(2*r+4))) %*% Matrix(psi2r6.ns, ncol=1), type="1")
if (!scv)
{
B1 <- 2*(2*pi)^(-d)*OF(2*r+4)*prod(d+2*(0:(r+2)))
B2 <- -(d+2*r+2)*2^(-d/2-r)*(2*pi)^(-d)*det(S)^(-1/2)*OF(2*r+4)*nu(r=r+2, A=Sinv)
}
else
{
B1 <- 2^(-2*r-3)*(4*pi)^(-d)*OF(2*r+6)*prod(d+2*(0:(r+2)))
B2 <- -(d+2*r+2)*2^(-2*r-4)*(4*pi)^(-d)*det(S)^(-1/2)*OF(2*r+4)*nu(r=r+2, A=Sinv)
B3 <- 4*B3
}
g2r4 <- (2*B1/(-B2 + sqrt(B2^2 + 4*B1*B3)))^(1/(d+2*r+8))*n^(-1/(d+2*r+8))
}
else if (nstage==2)
{
psi2r8.ns <- psins(r=2*r+8, Sigma=S)
B3 <- norm(Matrix(psi2r8.ns, nrow=1) %*% (Matrix(vec(diag(d)), ncol=1) %*% Matrix(vec(diag(d)), nrow=1) %x% Diagonal(d^(2*r+6))) %*% Matrix(psi2r8.ns, ncol=1), type="1")
if (!scv)
{
B1 <- 2*(2*pi)^(-d)*OF(2*r+6)*prod(d+2*(0:(r+3)))
B2 <- -(d+2*r+4)*2^(-d/2-r+2)*(2*pi)^(-d)*det(S)^(-1/2)*OF(2*r+6)*nu(r=r+4, A=Sinv)
}
else
{
B1 <- 2^(-2*r-5)*(4*pi)^(-d)*OF(2*r+6)*prod(d+2*(0:(r+3)))
B2 <- -(d+2*r+4)*2^(-2*r-6)*(4*pi)^(-d)*det(S)^(-1/2)*OF(2*r+6)*nu(r=r+4, A=Sinv)
B3 <- 4*B3
}
g2r6.ns <- (2*B1/(-B2 + sqrt(B2^2 + 4*B1*B3)))^(1/(d+2*r+8))*n^(-1/(d+2*r+8))
L0 <- dmvnorm.mixt(x=rep(0,d), mus=rep(0,d), Sigmas=diag(d), props=1)
if (binned)
eta2r6 <- drop(kfe(x=x, G=g2r6.ns^2*diag(d), inc=1, binned=binned, deriv.order=2*r+6, add.index=FALSE, verbose=verbose) %*% vec(diag(d^(r+3))))
else
eta2r6 <- Qr(x=x, deriv.order=2*r+6, Sigma=g2r6.ns^2*diag(d), inc=1)
A1 <- cf[1]*(2*d+4*r+8)*L0^2*OF(2*r+4)*nu(r=r+2, A=diag(d))
A2 <- cf[2]*(d+2*r+2)*L0*OF(2*r+4)*eta2r6
A3 <- cf[3]*eta2r6^2
g2r4 <- (2*A1/((-A2+ sqrt(A2^2 +4*A1*A3))*n))^(1/(d+2*r+6))
}
return(g2r4)
}
##############################################################################
## Unconstrained pilot selector for derivatives r>0 from Chacon & Duong (2011)
## Generalisation of Gunconstr for r>0
##############################################################################
Gdunconstr <- function(x, d, r, n, nstage=1, verbose, binned=FALSE, scv=FALSE, optim.fun="optim")
{
if (scv) cf <- c(2^(-(d/2+r+2)), 2)
else cf <- c(1,1)
S <- var(x)
S12 <- matrix.sqrt(S)
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
if (nstage==2)
{
# G2r6.NR <- Gns(r=2*r+6,n=n,Sigma=S, scv=scv)
# vecPsi2r6 <- kfe(x=x, G=G2r6.NR, binned=binned, deriv.order=2*r+6, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
x.star <- pre.sphere(x)
G2r6.NR <- Gns(r=2*r+6,n=n,Sigma=var(x.star), scv=scv)
vecPsi2r6 <- kfe(x=x.star, G=G2r6.NR, binned=binned, deriv.order=2*r+6, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
dls <- (0:(d^2-1))*d^(2*r+4)
AB2 <- function(vechG)
{
G <- invvech(vechG) %*% invvech(vechG)
Ginv <- chol2inv(chol(G))
## direct computation
v1 <- n^(-1)*det(G)^(-1/2)*(2*pi)^(-d/2)*(-1)^(r+2)*OF(2*r+4)* Sdrv(d=d,r=2*r+4,v=Kpow(vec(Ginv),r+2))
v2 <- numeric(d^(2*r+4))
for(k in 1:d^(2*r+4)) { v2[k] <- sum(vec(G)*vecPsi2r6[dls+k]) }
AB <- cf[1]*v1 + cf[2]*(1/2)*v2
AB2.val <- sum(AB^2)
return(AB2.val)
}
Gstart <- Gns(r=2*r+4, n=n, Sigma=var(x.star), scv=scv)
Gstart <- matrix.sqrt(Gstart)
if (optim.fun=="nlm")
{
result <- nlm(p=vech(Gstart), f=AB2, print.level=2*as.logical(verbose))
G2r4 <- result$estimate
}
else if (optim.fun=="optim")
{
result <- optim(par=vech(Gstart), fn=AB2, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
G2r4 <- result$par
}
G2r4 <- invvech(G2r4) %*% invvech(G2r4)
G2r4 <- S12 %*% G2r4 %*% S12
}
else
{
G2r4 <- Gns(r=2*r+4, n=n, Sigma=S, scv=scv)
}
return(G2r4)
}
##############################################################################
## Estimate psi functionals using 1-stage plug-in
##
## Parameters
## x.star - pre-transformed data points
## pilot - "amse" = different AMSE pilot bandwidths
## - "samse" = optimal SAMSE pilot bandwidth
##
## Returns
## estimated psi functionals
###############################################################################
psifun1 <- function(x.star, pilot="samse", binned, bin.par, deriv.order=0, verbose=FALSE)
{
d <- ncol(x.star)
r <- deriv.order
S.star <- var(x.star)
n <- nrow(x.star)
## pilots are based on (2r+4)-th order derivatives
## compute 1 pilot for SAMSE
if (pilot=="samse")
{
g.star <- gsamse(S.star, n, 4)
psihat.star <- kfe(x=x.star, G=g.star^2*diag(d), deriv.order=4, deriv.vec=TRUE, binned=binned, add.index=TRUE, verbose=verbose)
}
## compute 5 different pilots for AMSE
else if ((pilot=="amse") & (d==2))
{
derivt4 <- dmvnorm.deriv(x=rep(0,d), deriv.order=4, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
derivt4.vec <- dmvnorm.deriv(x=rep(0,d), deriv.order=4, add.index=TRUE, deriv.vec=TRUE, only.index=TRUE)
RK31 <- 15/(64*pi)
psi00 <- psins(r=c(0,0), Sigma=S.star)
psihat.star <- vector()
g.star <- vector()
for (k in 1:nrow(derivt4))
{
r <- derivt4[k,]
psi1 <- psins(r=r + 2*elem(1, 2), Sigma=S.star)
psi2 <- psins(r=r + 2*elem(2, 2), Sigma=S.star)
## odd order
if (prod(r) == 3) g.star[k] <- gamse.odd.2d(r=4, n, psi1, psi2, psi00, RK31)
## even order
else g.star[k] <- gamse.even.2d(r=4, n, psi1, psi2)[k]
psihat.star[k] <- kfe.scalar(x=x.star, deriv.order=r, g=g.star[k], binned=binned, bin.par=bin.par)
}
## create replicated form of psihat
psihat.star.vec <- rep(0, nrow(derivt4.vec))
for (k in 1:nrow(derivt4.vec))
psihat.star.vec[k] <- psihat.star[which.mat(r=derivt4.vec[k,], mat=derivt4)]
psihat.star <- list(psir=psihat.star.vec, deriv.ind=derivt4.vec)
}
return(psihat.star)
}
###############################################################################
## Estimate psi functionals using 2-stage plug-in
##
## Parameters
## x - pre-transformed data points
## pilot - "amse" - different AMSE pilot
## - "samse" - SAMSE pilot
## Returns
## estimated psi functionals
###############################################################################
psifun2 <- function(x.star, pilot="samse", binned, bin.par, deriv.order=0, verbose=FALSE)
{
d <- ncol(x.star)
r <- deriv.order
S.star <- var(x.star)
n <- nrow(x.star)
## pilots are based on (2r+4)-th order derivatives
## compute 1 pilot for SAMSE
if (pilot=="samse")
{
g6.star <- gsamse(S.star, n=n, modr=6)
psihat6.star <- kfe(x=x.star, G=g6.star^2*diag(d), deriv.order=6, deriv.vec=TRUE, binned=binned, add.index=FALSE, verbose=verbose)
g.star <- gsamse(S.star, n=n, modr=4, nstage=2, psihat=psihat6.star)
psihat.star <- kfe(x=x.star, G=g.star^2*diag(d), deriv.order=4, deriv.vec=TRUE, binned=binned, add.index=TRUE, verbose=verbose)
}
## compute different pilots for AMSE
else if ((pilot=="amse") & (d==2))
{
derivt4 <- dmvnorm.deriv(x=rep(0,d), deriv.order=4, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
derivt4.vec <- dmvnorm.deriv(x=rep(0,d), deriv.order=4, add.index=TRUE, deriv.vec=TRUE, only.index=TRUE)
derivt6 <- dmvnorm.deriv(x=rep(0,d), deriv.order=6, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
RK31 <- 15/(64*pi)
RK51 <- 945/(256*pi)
RK33 <- 225/(256*pi)
psi00 <- psins(r=rep(0,d), Sigma=S.star)
psihat6.star <- vector()
g6.star <- vector()
psihat.star <- vector()
g.star <- vector()
for (k in 1:nrow(derivt6))
{
r <- derivt6[k,]
psi1 <- psins(r=r + 2*elem(1, 2), Sigma=S.star)
psi2 <- psins(r=r + 2*elem(2, 2), Sigma=S.star)
if (prod(r) == 5)
g6.star[k] <- gamse.odd.2d(r=6, n, psi1, psi2, psi00, RK51)
else if (prod(r) == 9)
g6.star[k] <- gamse.odd.2d(r=6, n, psi1, psi2, psi00, RK33)
else
g6.star[k] <- gamse.even.2d(r=6, n, psi1, psi2)[k]
psihat6.star[k] <- kfe.scalar(x=x.star, deriv.order=r, g=g6.star[k], binned=binned, bin.par=bin.par)
}
## pilots are based on 4th order derivatives using 6th order psi functionals
## computed above 'psihat6.star'
for (k in 1:nrow(derivt4))
{
r <- derivt4[k,]
psi1 <- psihat6.star[7 - (r + 2*elem(1,2))[1]]
psi2 <- psihat6.star[7 - (r + 2*elem(2,2))[1]]
if (prod(r) == 3)
g.star[k] <- gamse.odd.2d(r=4, n, psi1, psi2, psi00, RK31)
else
g.star[k] <- gamse.even.2d(r=4, n, psi1, psi2)[k]
psihat.star[k] <- kfe.scalar(x=x.star, deriv.order=r, g=g.star[k], binned=binned, bin.par=bin.par)
}
## create replicated form of psihat
psihat.star.vec <- rep(0, nrow(derivt4.vec))
for (k in 1:nrow(derivt4.vec))
psihat.star.vec[k] <- psihat.star[which.mat(r=derivt4.vec[k,], mat=derivt4)]
psihat.star <- list(psir=psihat.star.vec, deriv.ind=derivt4.vec)
}
return(psihat.star)
}
#############################################################################
## Estimate psi functionals for 6-variate data using 1-stage plug-in
## with unconstrained pilot
##
## Parameters
## x - data points
## Sd4, Sd6 - symmetrizer matrices of order 4 and 6
##
## Returns
## estimated psi functionals
#############################################################################
psifun1.unconstr <- function(x, binned, bgridsize, deriv.order=0, verbose=FALSE)
{
n <- nrow(x)
r <- deriv.order
S <- var(x)
## stage 1 of plug-in
G2r4 <- Gns(r=2*r+4,n=n,Sigma=S)
vecPsi2r4 <- kfe(x=x, G=G2r4, deriv.order=2*r+4, binned=binned, bgridsize=bgridsize, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
return (vecPsi2r4)
}
#############################################################################
## Estimate psi functionals for 6-variate data using 2-stage plug-in
## with unconstrained pilot
##
## Parameters
## x - data points
## Sd4, Sd6 - symmetrizer matrices of order 4 and 6
##
## Returns
## estimated psi functionals
############################################################################
psifun2.unconstr <- function(x, binned, bgridsize, deriv.order=0, verbose=FALSE, optim.fun="optim")
{
d <- ncol(x)
n <- nrow(x)
S <- var(x)
S12 <- matrix.sqrt(S)
r <- deriv.order
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
## stage 1 of plug-in
x.star <- pre.sphere(x)
G2r6 <- Gns(r=2*r+6,n=n,Sigma=var(x.star))
vecPsi2r6 <- kfe(x=x.star, G=G2r6, binned=binned, bgridsize=bgridsize, deriv.order=2*r+6, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
## asymptotic squared bias for r = 4 for MSE-optimal G
D2r4phi0 <- DrL0(d=d, r=2*r+4)
Id2r4 <- diag(d^(2*r+4))
AB2 <- function(vechG)
{
rr <- 2*r+4
G <- invvech(vechG)%*%invvech(vechG)
G12 <- matrix.sqrt(G)
Ginv12 <- chol2inv(chol(G12))
AB <- n^(-1)*det(Ginv12)*(Kpow(A=Ginv12,pow=rr)%*%D2r4phi0)+(1/2)*(t(vec(G))%x%Id2r4) %*% vecPsi2r6
AB2.val <- sum(AB^2)
return (AB2.val)
}
Gstart <- Gns(r=2*r+4,n=n,Sigma=var(x.star))
Gstart <- matrix.sqrt(Gstart)
if (optim.fun=="nlm")
{
res <- nlm(p=vech(Gstart), f=AB2, print.level=2*as.logical(verbose))
G2r4 <- res$estimate
}
else if (optim.fun=="optim")
{
res <- optim(vech(Gstart), AB2, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
G2r4 <- res$par
}
G2r4 <- invvech(G2r4)%*%invvech(G2r4)
G2r4 <- S12 %*% G2r4 %*% S12
## stage 2 of plug-in
vecPsi2r4 <- kfe(x=x, G=G2r4, binned=binned, bgridsize=bgridsize, deriv.order=2*r+4, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
return (vecPsi2r4)
}
#############################################################################
## Plug-in bandwidth selectors
#############################################################################
############################################################################
## Computes plug-in full bandwidth matrix - 2 to 6 dim
##
## Parameters
## x - data points
## Hstart - initial value for minimisation
## nstage - number of plug-in stages (1 or 2)
## pilot - "amse" - different AMSE pilot
## - "samse" - SAMSE pilot
## - "unconstr" - unconstrained pilot
## pre - "scale" - pre-scaled data
## - "sphere"- pre-sphered data
##
## Returns
## Plug-in full bandwidth matrix
###############################################################################
hpi <- function(x, nstage=2, binned=TRUE, bgridsize, deriv.order=0)
{
## 1-d selector is taken from KernSmooth's dpik
d <- 1
if (missing(bgridsize)) bgridsize <- default.gridsize(d)
if (deriv.order==0) h <- dpik(x=x, level=nstage, gridsize=bgridsize)
else
{
n <- length(x)
d <- 1
r <- deriv.order
K2r4 <- dnorm.deriv(x=0, mu=0, sigma=1, deriv.order=2*r+4)
K2r6 <- dnorm.deriv(x=0, mu=0, sigma=1, deriv.order=2*r+6)
m2 <- 1
mr <- psins.1d(r=2*r, sigma=1)
## formula for bias annihilating bandwidths from Wand & Jones (1995, p.70)
if (nstage==2)
{
psi2r8.hat <- psins.1d(r=2*r+8, sigma=sd(x))
gamse2r6 <- (2*K2r6/(-m2*psi2r8.hat*n))^(1/(2*r+9))
psi2r6.hat <- kfe.1d(x=x, g=gamse2r6, deriv.order=2*r+6, inc=1, binned=binned)
gamse2r4 <- (2*K2r4/(-m2*psi2r6.hat*n))^(1/(2*r+7))
psi2r4.hat <- kfe.1d(x=x, g=gamse2r4, deriv.order=2*r+4, inc=1, binned=binned)
}
else
{
psi2r6.hat <- psins.1d(r=2*r+6, sigma=sd(x))
gamse2r4 <- (2*K2r4/(-m2*psi2r6.hat*n))^(1/(2*r+7))
psi2r4.hat <- kfe.1d(x=x, g=gamse2r4, deriv.order=2*r+4, inc=1, binned=binned)
}
## formula form Wand & Jones (1995, p.49)
h <- ((2*r+1)*mr/(m2^2*psi2r4.hat*n))^(1/(2*r+5))
}
return(h)
}
Hpi <- function(x, nstage=2, pilot, pre="sphere", Hstart, binned, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim")
{
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d > 4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
if (!is.matrix(x)) x <- as.matrix(x)
if (missing(pilot)) { if (d==2 & r==0) pilot <- "samse" else pilot <- "dscalar" }
pilot1 <- match.arg(pilot, c("amse", "samse", "unconstr", "dunconstr", "dscalar"))
pre1 <- match.arg(pre, c("scale", "sphere"))
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
if (pilot1=="amse" & (d>2 | r>0)) stop("amse pilot selectors not defined for d>2 and/or r>0")
if ((pilot1=="samse" | pilot1=="unconstr") & r>0) stop("dscalar or dunconstr pilot selectors are better for derivatives r>0")
if (pilot1=="unconstr" & d>=6) stop("Unconstrained pilots are not implemented for d>6")
if (pre1=="scale")
{
x.star <- pre.scale(x)
S12 <- diag(sqrt(diag(var(x))))
Sinv12 <- chol2inv(chol(S12))
}
else if (pre1=="sphere")
{
x.star <- pre.sphere(x)
S12 <- matrix.sqrt(var(x))
Sinv12 <- chol2inv(chol(S12))
}
Idr <- diag(d^r)
RKr <- nu(r=r, diag(d))*2^(-d-r)*pi^(-d/2)
if (binned)
{
H.max <- (((d+8)^((d+6)/2)*pi^(d/2)*RKr)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x)
bin.par <- binning(x=x, bgridsize=bgridsize, H=H.max)
H.max.star <- (((d+8)^((d+6)/2)*pi^(d/2)*RKr)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
bin.par.star <- binning(x=x.star, bgridsize=bgridsize, H=H.max.star)
}
if (pilot1=="unconstr")
{
## psi4.mat is on data scale
if (nstage==1)
psi.fun <- (-1)^r*psifun1.unconstr(x=x, binned=binned, bgridsize=bgridsize, deriv.order=r, verbose=verbose)
else if (nstage==2)
psi.fun <- psifun2.unconstr(x=x, binned=binned, bgridsize=bgridsize, deriv.order=r, verbose=verbose)
psi2r4.mat <- (-1)^r*invvec(psi.fun)
## use normal reference bandwidth as initial condition
if (missing(Hstart)) Hstart <- Hns(x=x, deriv.order=r)
}
else if (pilot1=="dunconstr")
{
## G2r4 is on data scale
G2r4 <- Gdunconstr(x=x, d=d, r=r, n=n, nstage=nstage, verbose=verbose, binned=binned, optim.fun=optim.fun)
vecPsi2r4 <- kfe(x=x, G=G2r4, binned=binned, deriv.order=2*r+4, deriv.vec=TRUE, add.index=FALSE, verbose=verbose, bin.par=bin.par)
if (missing(Hstart)) Hstart <- Hns(x=x, deriv.order=r)
}
else if (pilot1=="dscalar")
{
## g2r4 is on pre-transformed data scale
g2r4 <- gdscalar(x=x.star, r=r, n=n, d=d, verbose=verbose, nstage=nstage, binned=binned)
G2r4 <- g2r4^2 * diag(d)
vecPsi2r4 <- kfe(x=x.star, G=G2r4, binned=binned, deriv.order=2*r+4, deriv.vec=TRUE, add.index=FALSE, verbose=verbose, bin.par=bin.par.star)
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
}
else
{
## psi4.mat is on pre-transformed data scale
if (nstage==1)
psi.fun <- psifun1(x.star, pilot=pilot, binned=binned, bin.par=bin.par.star, deriv.order=r, verbose=verbose)$psir
else if (nstage==2)
psi.fun <- psifun2(x.star, pilot=pilot, binned=binned, bin.par=bin.par.star, deriv.order=r, verbose=verbose)$psir
psi2r4.mat <- invvec(psi.fun)
## use normal reference bandwidth as initial condition
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
else Hstart <- Sinv12 %*% Hstart %*% Sinv12
}
## PI is estimate of AMISE
pi.temp <- function(vechH)
{
H <- invvech(vechH) %*% invvech(vechH)
Hinv <- chol2inv(chol(H))
IdrvH <- Idr%x%vec(H)
int.var <- 1/(det(H)^(1/2)*n)*nur(r=r, A=Hinv, mu=rep(0,d), Sigma=diag(d))*2^(-d-r)*pi^(-d/2)
if (pilot1=="dunconstr" | pilot1=="dscalar")
pi.val <- int.var + (-1)^r*1/4*vecPsi2r4 %*% (vec(diag(d^r) %x% vec(H) %x% vec(H)))
else
pi.val <- int.var + (-1)^r*1/4* sum(diag(t(IdrvH) %*% psi2r4.mat %*% IdrvH))
pi.val <- drop(pi.val)
return(pi.val)
}
Hstart <- matrix.sqrt(Hstart)
if (optim.fun=="nlm")
{
result <- nlm(p=vech(Hstart), f=pi.temp, print.level=2*as.numeric(verbose))
H <- invvech(result$estimate) %*% invvech(result$estimate)
amise.star <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(vech(Hstart), pi.temp, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
H <- invvech(result$par) %*% invvech(result$par)
amise.star <- result$value
}
if (!(pilot1 %in% c("dunconstr","unconstr"))) H <- S12 %*% H %*% S12 ## back-transform
if (!amise) return(H)
else return(list(H = H, PI.star=amise.star))
}
###############################################################################
## Computes plug-in diagonal bandwidth matrix for 2 to 6-dim
##
## Parameters
## x - data points
## nstage - number of plug-in stages (1 or 2)
## pre - "scale" - pre-scaled data
## - "sphere"- pre-sphered data
##
## Returns
## Plug-in diagonal bandwidth matrix
###############################################################################
Hpi.diag <- function(x, nstage=2, pilot, pre="scale", Hstart, binned, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim")
{
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
RK <- (4*pi)^(-d/2)
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d > 4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
if (!is.matrix(x)) x <- as.matrix(x)
if (missing(pilot)) { if (d==2 & r==0) pilot <- "samse" else pilot <- "dscalar" }
pilot1 <- match.arg(pilot, c("amse", "samse", "unconstr", "dunconstr", "dscalar"))
pre1 <- match.arg(pre, c("scale", "sphere"))
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
if (pre1=="sphere") stop("Using pre-sphering won't give diagonal bandwidth matrix")
if (pilot1=="amse" & (d>2 | r>0)) stop("samse pilot selectors are better for higher dimensions and/or deriv.order>0")
if (pilot1=="samse" & r>0) stop("dscalar or dunconstr pilot selectors are better for derivatives r>0")
if (pilot1=="unconstr" | pilot1=="dunconstr") stop("Unconstrained pilot selectors are not suitable for Hpi.diag")
if (pre1=="scale")
{
x.star <- pre.scale(x)
S12 <- diag(sqrt(diag(var(x))))
Sinv12 <- chol2inv(chol(S12))
}
else if (pre1=="sphere")
{
x.star <- pre.sphere(x)
S12 <- matrix.sqrt(var(x))
Sinv12 <- chol2inv(chol(S12))
}
if (binned)
{
H.max <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
bin.par <- binning(x=x.star, bgridsize=bgridsize, H=H.max)
}
Idr <- diag(d^r)
if (pilot1=="amse" | pilot1=="samse")
{
if (nstage==1)
psi.fun <- psifun1(x.star, pilot=pilot, binned=binned, bin.par=bin.par, deriv.order=r, verbose=verbose)$psir
else if (nstage==2)
psi.fun <- psifun2(x.star, pilot=pilot, binned=binned, bin.par=bin.par, deriv.order=r, verbose=verbose)$psir
psi2r4.mat <- invvec(psi.fun)
}
else if (pilot1=="dscalar")
{
g2r4 <- gdscalar(x=x.star, r=r, n=n, d=d, verbose=verbose, nstage=nstage, binned=binned)
G2r4 <- g2r4^2 * diag(d)
vecPsi2r4 <- kfe(x=x.star, G=G2r4, binned=binned, deriv.order=2*r+4, deriv.vec=TRUE, add.index=FALSE, verbose=verbose)
}
if (d==2 & r==0 & (pilot1=="amse" | pilot1=="samse"))
{
## diagonal bandwidth matrix for 2-dim has exact formula
psi40 <- psi.fun[1]
psi22 <- psi.fun[6]
psi04 <- psi.fun[16]
s1 <- sd(x[,1])
s2 <- sd(x[,2])
h1 <- (psi04^(3/4)*RK/(psi40^(3/4)*(sqrt(psi40*psi04)+psi22)*n))^(1/6)
h2 <- (psi40/psi04)^(1/4) * h1
H <- diag(c(s1^2*h1^2, s2^2*h2^2))
psimat4.D <- invvech(c(psi40, psi22, psi04))
amise.star <- drop(n^(-1)*RK*(h1*h2)^(-1) + 1/4*c(h1,h2)^2 %*% psimat4.D %*% c(h1,h2)^2)
}
else
{
## PI is estimate of AMISE
pi.temp <- function(diagH)
{
H <- diag(diagH) %*% diag(diagH)
Hinv <- chol2inv(chol(H))
IdrvH <- Idr%x%vec(H)
int.var <- 1/(det(H)^(1/2)*n)*nu(r=r, Hinv)*2^(-d-r)*pi^(-d/2)
if (pilot1=="dscalar")
pi.val <- int.var + (-1)^r*1/4*vecPsi2r4 %*% (vec(diag(d^r) %x% vec(H) %x% vec(H)))
else
pi.val <- int.var + (-1)^r*1/4* sum(diag(t(IdrvH) %*% psi2r4.mat %*% IdrvH))
return(drop(pi.val))
}
## use normal reference bandwidth as initial condition
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
else Hstart <- Sinv12 %*% Hstart %*% Sinv12
Hstart <- matrix.sqrt(Hstart)
if (optim.fun=="nlm")
{
result <- nlm(p=diag(Hstart), f=pi.temp, print.level=2*as.numeric(verbose))
H <- diag(result$estimate^2)
amise.star <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(diag(Hstart), pi.temp, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
H <- diag(result$par) %*% diag(result$par)
amise.star <- result$value
}
H <- S12 %*% H %*% S12 ## back-transform
}
if (!amise) return(H)
else return(list(H = H, PI.star=amise.star))
}
###############################################################################
## Cross-validation bandwidth selectors
###############################################################################
###############################################################################
## Computes the least squares cross validation LSCV function for 2 to 6 dim
##
## Parameters
## x - data values
## H - bandwidth matrix
##
## Returns
## LSCV(H)
###############################################################################
lscv.1d <- function(x, h, binned, bin.par, deriv.order=0)
{
r <- deriv.order
lscv1 <- kfe.1d(x=x, g=sqrt(2)*h, inc=1, binned=binned, bin.par=bin.par, deriv.order=2*r)
lscv2 <- kfe.1d(x=x, g=h, inc=0, binned=binned, bin.par=bin.par, deriv.order=2*r)
return((-1)^r*(lscv1 - 2*lscv2))
}
lscv.mat <- function(x, H, binned=FALSE, bin.par, bgridsize, deriv.order=0)
{
r <- deriv.order
d <- ncol(x)
n <- nrow(x)
if (!binned)
{
lscv1 <- Qr(x=x, deriv.order=2*r, Sigma=2*H, inc=1)
lscv2 <- Qr(x=x, deriv.order=2*r, Sigma=H, inc=0)
lscv <- drop(lscv1 - 2*lscv2)
}
else
{
lscv1 <- kfe(x=x, G=2*H, inc=1, binned=binned, bin.par=bin.par, bgridsize=bgridsize, deriv.order=2*r, add.index=FALSE)
lscv2 <- kfe(x=x, G=H, inc=0, binned=binned, bin.par=bin.par, bgridsize=bgridsize, deriv.order=2*r, add.index=FALSE)
lscv <- (-1)^2*drop((lscv1 - 2*lscv2) %*% vec(diag(d^r)))
}
return(lscv)
}
###############################################################################
## Finds the bandwidth matrix that minimises LSCV for 2 to 6 dim
##
## Parameters
## x - data values
## Hstart - initial bandwidth matrix
##
## Returns
## H_LSCV
###############################################################################
hlscv <- function(x, binned=TRUE, bgridsize, amise=FALSE, deriv.order=0, bw.ucv=TRUE)
{
## adapted from versions supplied by J.E. Chacon 19/02/2021
if (any(duplicated(x)))
warning("Data contain duplicated values: LSCV is not well-behaved in this case")
n <- length(x)
d <- 1
r <- deriv.order
hnorm <- sqrt((4/(n*(d + 2)))^(2/(d + 4)) * var(x))
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (missing(bgridsize)) if (bw.ucv) bgridsize <- 10000 else bgridsize <- 10001
## use stats::bw.ucv function
if (bw.ucv)
{
difs <- dist(x)
lower <- min(difs[difs>0])
opt <- list()
opt$minimum <- stats::bw.ucv(x=x, nb=bgridsize, lower=lower)
opt$objective <- NA
}
else
{
difs <- x%*%t(rep(1,n))-rep(1,n)%*%t(x)
difs <- difs[lower.tri(difs)]
if (binned)
{
bin.par <- binning(x, bgridsize=bgridsize, h=hnorm)
lscv.1d.temp <- function(h) { return(lscv.1d(x=x, h=h, binned=binned, bin.par=bin.par, deriv.order=r)) }
}
else
{
if (r>0) stop("Unbinned hlscv not yet implemented for deriv.order>0")
edifs <- exp(-difs^2/2)
RK <- 1/(2*sqrt(pi))
lscv.1d.temp <- function(h)
{
lscv1 <- (1-1/n)*sum(edifs^(1/(2*h^2)))/(h*sqrt(2)*sqrt(2*pi))
lscv2 <- 2*sum(edifs^(1/h^2))/(h*sqrt(2*pi))
return(RK/(n*h)+2*(lscv1-lscv2)/(n^2-n))
}
}
lower <- min(difs[difs>0])
opt <- optimise(lscv.1d.temp, interval=c(lower, 2*hnorm), tol=.Machine$double.eps)
}
if (!amise) return(opt$minimum)
else return(list(h=opt$minimum, LSCV=opt$objective))
}
Hlscv <- function(x, Hstart, binned=FALSE, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim", trunc)
{
if (any(duplicated(x))) warning("Data contain duplicated values: LSCV is not well-behaved in this case")
if (!is.matrix(x)) x <- as.matrix(x)
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d>4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
## use normal reference selector as initial condn
Hnorm <- Hns(x=x, deriv.order=r)
if (missing(Hstart)) Hstart <- Hnorm
Hstart <- matrix.sqrt(Hstart)
if (binned) bin.par <- binning(x=x, H=Hnorm, bgridsize=bgridsize)
if (missing(trunc)) { if (deriv.order==0) trunc <- 1e10 else trunc <- 4 }
lscv.init <- lscv.mat(x=x, H=Hnorm, binned=binned, bin.par=bin.par, deriv.order=r)
lscv.mat.temp <- function(vechH)
{
H <- invvech(vechH) %*% invvech(vechH)
lscv <- lscv.mat(x=x, H=H, binned=binned, bin.par=bin.par, deriv.order=r)
if (det(H) < 1/trunc*det(Hnorm) | det(H) > trunc*det(Hnorm) | abs(lscv) > trunc*abs(lscv.init)) lscv <- lscv.init
return(lscv)
}
if (optim.fun=="nlm")
{
result <- nlm(p=vech(Hstart), f=lscv.mat.temp, print.level=2*as.numeric(verbose))
H <- invvech(result$estimate) %*% invvech(result$estimate)
amise.opt <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(vech(Hstart), lscv.mat.temp, method="Nelder-Mead", control=list(trace=as.numeric(verbose), REPORT=1))
H <- invvech(result$par) %*% invvech(result$par)
amise.opt <- result$value
}
if (!amise) return(H)
else return(list(H=H, LSCV=amise.opt))
}
###############################################################################
## Finds the diagonal bandwidth matrix that minimises LSCV for 2 to 6 dim
##
## Parameters
## x - data values
## Hstart - initial bandwidth matrix
##
## Returns
## H_LSCV,diag
###############################################################################
Hlscv.diag <- function(x, Hstart, binned=FALSE, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim", trunc)
{
if (any(duplicated(x))) warning("Data contain duplicated values: LSCV is not well-behaved in this case")
if (!is.matrix(x)) x <- as.matrix(x)
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d>4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
Hnorm <- Hns(x=x, deriv.order=r)
if (missing(Hstart)) Hstart <- Hnorm
## don't truncate optimisation for deriv.order==0
if (missing(trunc)) { if (deriv.order==0) trunc <- 1e10 else trunc <- 4 }
## linear binning
if (binned)
{
bin.par <- binning(x=x, bgridsize=bgridsize, H=Hnorm)
}
lscv.init <- lscv.mat(x=x, H=Hnorm, binned=binned, bin.par=bin.par, deriv.order=r)
lscv.mat.temp <- function(diagH)
{
H <- diag(diagH^2)
lscv <- lscv.mat(x=x, H=H, binned=binned, bin.par=bin.par, deriv.order=r)
if (det(H) < 1/trunc*det(Hnorm) | det(H) > trunc*det(Hnorm) | abs(lscv) > trunc*abs(lscv.init)) lscv <- lscv.init
return(lscv)
}
Hstart <- matrix.sqrt(Hstart)
if (optim.fun=="nlm")
{
result <- nlm(p=diag(Hstart), f=lscv.mat.temp, print.level=2*as.numeric(verbose))
H <- diag(result$estimate^2)
amise.opt <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(diag(Hstart), lscv.mat.temp, method="Nelder-Mead", control=list(trace=as.numeric(verbose), REPORT=1))
H <- diag(result$par^2)
amise.opt <- result$value
}
if (!amise) return(H)
else return(list(H=H, LSCV=amise.opt))
}
## aliases using "UCV" instead of "LSCV"
hucv <- function(...) { hlscv(...) }
Hucv <- function(...) { Hlscv(...) }
Hucv.diag <- function(...) { Hlscv.diag(...) }
###############################################################################
## Computes the biased cross validation BCV function for 2-dim
##
## Parameters
## x - data values
## H1, H2 - bandwidth matrices
##
## Returns
## BCV(H)
###############################################################################
bcv.mat <- function(x, H1, H2, binned=FALSE)
{
n <- nrow(x)
d <- 2
psi <- kfe(x, G=H2, deriv.order=4, add.index=TRUE, deriv.vec=TRUE, inc=0, binned=binned)
psi40 <- psi$psir[1]
psi31 <- psi$psir[2]
psi22 <- psi$psir[4]
psi13 <- psi$psir[8]
psi04 <- psi$psir[16]
coeff <- c(1, 2, 1, 2, 4, 2, 1, 2, 1)
psi.fun <- c(psi40, psi31, psi22, psi31, psi22, psi13, psi22, psi13,psi04)/(n*(n-1))
psi4.mat <- matrix(coeff * psi.fun, ncol=3, nrow=3)
RK <- (4*pi)^(-d/2)
bcv <- drop(n^(-1)*det(H1)^(-1/2)*RK + 1/4*t(vech(H1)) %*% psi4.mat %*% vech(H1))
return(list(bcv=bcv, psimat=psi4.mat))
}
###############################################################################
## Find the bandwidth matrix that minimises the BCV for 2-dim
##
## Parameters
## x - data values
## whichbcv - 1 = BCV1
## - 2 = BCV2
## Hstart - initial bandwidth matrix
##
## Returns
## H_BCV
###############################################################################
Hbcv <- function(x, whichbcv=1, Hstart, binned=FALSE, amise=FALSE, verbose=FALSE)
{
n <- nrow(x)
d <- ncol(x)
RK <- (4*pi)^(-d/2)
if(!is.matrix(x)) x <- as.matrix(x)
## use normal reference b/w matrix for bounds
k <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*n*gamma((d+8)/2)*(d+2)))^(2/(d+4))
Hmax <- k * abs(var(x))
up.bound <- Hmax
if (missing(Hstart)) Hstart <- matrix.sqrt(0.9*Hmax)
bcv.mat.temp <- function(vechH)
{
H <- invvech(vechH) %*% invvech(vechH)
bcv <- bcv.mat(x=x, H1=H, H2=whichbcv*H, binned=binned)$bcv
return(bcv)
}
result <- optim(vech(Hstart), bcv.mat.temp, method="L-BFGS-B", upper=vech(matrix.sqrt(up.bound)), lower=-vech(matrix.sqrt(up.bound)), control=list(trace=as.numeric(verbose), REPORT=1))
H <- invvech(result$par) %*% invvech(result$par)
amise.opt <- result$value
if (!amise) return(H)
else return(list(H = H, BCV=amise.opt))
}
###############################################################################
## Find the diagonal bandwidth matrix that minimises the BCV for 2-dim
##
## Parameters
## x - data values
## whichbcv - 1 = BCV1
## - 2 = BCV2
## Hstart - initial bandwidth matrix
##
## Returns
## H_BCV, diag
###############################################################################
Hbcv.diag <- function(x, whichbcv=1, Hstart, binned=FALSE, amise=FALSE, verbose=FALSE)
{
n <- nrow(x)
d <- ncol(x)
RK <- (4*pi)^(-d/2)
if(!is.matrix(x)) x <- as.matrix(x)
## use maximally smoothed b/w matrix for bounds
k <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*n*gamma((d+8)/2)*(d+2)))^(2/(d+4))
Hmax <- k * abs(var(x))
up.bound <- diag(Hmax)
if (missing(Hstart)) Hstart <- 0.9*matrix.sqrt(Hmax)
bcv.mat.temp <- function(diagH)
{
H <- diag(diagH) %*% diag(diagH)
return(bcv.mat(x, H, whichbcv*H, binned=binned)$bcv)
}
result <- optim(diag(Hstart), bcv.mat.temp, method="L-BFGS-B", upper=sqrt(up.bound), control=list(trace=as.numeric(verbose), REPORT=1))
H <- diag(result$par) %*% diag(result$par)
amise.opt <- result$value
if (!amise) return(H)
else return(list(H = H, BCV=amise.opt))
}
###############################################################################
## Estimate scalar g_AMSE pilot bandwidth for SCV for 2 to 6 dim
##
## Parameters
## Sigma.star - scaled/ sphered variance matrix
## Hamise - (estimate) of H_AMISE
## n - sample size
##
## Returns
## g_AMSE pilot bandwidth
###############################################################################
Theta6.elem <- function(d)
{
Theta6.mat <- list()
Theta6.mat[[d]] <- list()
for (i in 1:d)
Theta6.mat[[i]] <- list()
for (i in 1:d)
for (j in 1:d)
{
temp <- numeric()
for (k in 1:d)
for (ell in 1:d)
temp <- rbind(temp, elem(i,d)+2*elem(k,d)+2*elem(ell,d)+elem(j,d))
Theta6.mat[[i]][[j]] <- temp
}
return(Theta6.mat)
}
gamse.scv <- function(x.star, d, Sigma.star, Hamise, n, binned=FALSE, bin.par, bgridsize, verbose=FALSE, nstage=1, Theta6=FALSE)
{
if (nstage==0)
{
psihat6.star <- psins(r=6, Sigma=Sigma.star, deriv.vec=TRUE)
}
else if (nstage==1)
{
g6.star <- gsamse(Sigma.star, n, 6)
G6.star <- g6.star^2 * diag(d)
if (Theta6) psihat6.star <- kfe(x=x.star, bin.par=bin.par, deriv.order=6, G=G6.star, deriv.vec=FALSE, add.index=FALSE, binned=binned, bgridsize=bgridsize, verbose=verbose)
else psihat6.star <- kfe(x=x.star, bin.par=bin.par, deriv.order=6, G=G6.star, deriv.vec=TRUE, add.index=FALSE, binned=binned, bgridsize=bgridsize, verbose=verbose)
}
if (Theta6)
{
derivt6 <- dmvnorm.deriv(x=rep(0,d), deriv.order=6, add.index=TRUE, deriv.vec=FALSE, only.index=TRUE)
Theta6.mat <- matrix(0, ncol=d, nrow=d)
Theta6.mat.ind <- Theta6.elem(d)
for (i in 1:d)
for (j in 1:d)
{
temp <- Theta6.mat.ind[[i]][[j]]
temp.sum <- 0
for (k in 1:nrow(temp))
temp.sum <- temp.sum + psihat6.star[which.mat(temp[k,], derivt6)]
Theta6.mat[i,j] <- temp.sum
}
eye3 <- diag(d)
D4 <- dupl(d)$d
trHamise <- tr(Hamise)
## required constants - see thesis
Cmu1 <- 1/2*t(D4) %*% vec(Theta6.mat %*% Hamise)
Cmu2 <- 1/8*(4*pi)^(-d/2) * (2*t(D4)%*% vec(Hamise) + trHamise * t(D4) %*% vec(eye3))
num <- 2 * (d+4) * sum(Cmu2*Cmu2)
den <- -(d+2) * sum(Cmu1*Cmu2) + sqrt((d+2)^2 * sum(Cmu1*Cmu2)^2 + 8*(d+4)*sum(Cmu1*Cmu1) * sum(Cmu2*Cmu2))
gamse <- (num/(den*n))^(1/(d+6))
}
else
{
## updated constants using Chacon & Duong (2010) notation
Cmu1Cmu1 <- drop(1/4*psihat6.star %*% (Hamise %x% diag(d^4) %x% Hamise) %*% psihat6.star)
Cmu1Cmu2 <- 3/4*(4*pi)^(-d/2)*drop(vec(Hamise %x% diag(d) %x% Hamise) %*% psihat6.star)
Cmu2Cmu2 <- 1/64*(4*pi)^(-d)*(4*tr(Hamise%*%Hamise) + (d+8)*tr(Hamise)^2)
num <- 2 * (d+4) * Cmu2Cmu2
den <- -(d+2) * Cmu1Cmu2 + sqrt((d+2)^2 * Cmu1Cmu2^2 + 8*(d+4)*Cmu1Cmu1 * Cmu2Cmu2)
gamse <- (num/(den*n))^(1/(d+6))
}
return(gamse)
}
###############################################################################
## Estimate unconstrained G_AMSE pilot bandwidth for SCV for 2 to 6 dim
## Code by J.E. Chacon
##
## Returns
## G_AMSE pilot bandwidth
###############################################################################
Gunconstr.scv <- function(x, binned=FALSE, bin.par, bgridsize, verbose=FALSE, nstage=1, optim.fun="optim")
{
d <- ncol(x)
n <- nrow(x)
S <- var(x)
S12 <- matrix.sqrt(S)
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
## stage 1 of plug-in
x.star <- pre.sphere(x)
if (nstage==1)
{
G6 <- Gns(r=6,n=n,Sigma=var(x.star), scv=TRUE)
psihat6 <- kfe(x=x.star, deriv.order=6, G=G6, deriv.vec=TRUE, add.index=FALSE, binned=binned, bgridsize=bgridsize, verbose=verbose)
}
else if (nstage==0)
{
psihat6 <- psins(r=6, Sigma=var(x.star), deriv.vec=TRUE)
}
## constants for normal reference
D4phi0 <- DrL0(d=d,r=4)
Id4 <- diag(d^4)
## asymptotic squared bias for r = 4
AB2 <- function(vechG)
{
G <- invvech(vechG)%*%invvech(vechG)
G12 <- matrix.sqrt(G)
Ginv12 <- chol2inv(chol(G12))
AB <- n^(-1)*det(Ginv12)*(Kpow(A=Ginv12,pow=4)%*%D4phi0)*2^(-(d+4)/2) + (t(vec(G))%x%Id4)%*%psihat6
AB2.val <- sum(AB^2)
return (AB2.val)
}
Gstart <- Gns(r=4,n=n,Sigma=S,scv=TRUE)
Gstart <- matrix.sqrt(Gstart)
if (optim.fun=="nlm")
{
result <- nlm(p=vech(Gstart), f=AB2, print.level=2*as.logical(verbose))
G4 <- result$estimate
}
else if (optim.fun=="optim")
{
result <- optim(vech(Gstart), AB2, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
G4 <- result$par
}
G4 <- invvech(G4)%*%invvech(G4)
G4 <- S12 %*% G4 %*% S12
return(G4)
}
###############################################################################
## Computes the smoothed cross validation function for 2 to 6 dim
##
## Parameters
## x - data values
## H - bandwidth matrix
## G - pilot bandwidth matrix
##
## Returns
## SCV(H)
###############################################################################
scv.1d <- function(x, h, g, binned=TRUE, bin.par, inc=1, deriv.order=0)
{
r <- deriv.order
if (!missing(x)) n <- length(x)
if (!missing(bin.par)) n <- sum(bin.par$counts)
scv1 <- kfe.1d(x=x, deriv.order=2*r, bin.par=bin.par, g=sqrt(2*h^2+2*g^2), binned=binned, inc=inc)
scv2 <- kfe.1d(x=x, deriv.order=2*r, bin.par=bin.par, g=sqrt(h^2+2*g^2), binned=binned, inc=inc)
scv3 <- kfe.1d(x=x, deriv.order=2*r, bin.par=bin.par, g=sqrt(2*g^2), binned=binned, inc=inc)
bias2 <- (-1)^r*(scv1 - 2*scv2 + scv3)
if (bias2 < 0) bias2 <- 0
scv <- (n*h)^(-1)*(4*pi)^(-1/2)*2^(-r)*OF(2*r) + bias2
return(scv)
}
scv.mat <- function(x, H, G, binned=FALSE, bin.par, bgridsize, verbose=FALSE, deriv.order=0)
{
d <- ncol(x)
n <- nrow(x)
r <- deriv.order
vId <- vec(diag(d))
Hinv <- chol2inv(chol(H))
if (!binned)
{
scv1 <- Qr(x=x, deriv.order=2*r, Sigma=2*H+2*G, inc=1)
scv2 <- Qr(x=x, deriv.order=2*r, Sigma=H+2*G, inc=1)
scv3 <- Qr(x=x, deriv.order=2*r, Sigma=2*G, inc=1)
bias2 <- (-1)^r*(scv1 - 2*scv2 + scv3)
if (bias2 < 0) bias2 <- 0
}
else
{
scv1 <- kfe(x=x, G=2*H + 2*G, deriv.order=2*r, inc=1, binned=binned, bin.par=bin.par, bgridsize=bgridsize, verbose=verbose, add.index=FALSE)
scv2 <- kfe(x=x, G=H + 2*G, deriv.order=2*r, inc=1, binned=binned, bin.par=bin.par, bgridsize=bgridsize, verbose=verbose, add.index=FALSE)
scv3 <- kfe(x=x, G=2*G, deriv.order=2*r, inc=1, binned=binned, bin.par=bin.par, bgridsize=bgridsize, verbose=verbose, add.index=FALSE)
bias2 <- drop((-1)^r*Kpow(vId,r) %*% (scv1 - 2*scv2 + scv3))
if (bias2 < 0) bias2 <- 0
}
scvmat <- 1/(det(H)^(1/2)*n)*nur(r=r, A=Hinv, mu=rep(0,d), Sigma=diag(d), type="direct")*2^(-d-r)*pi^(-d/2) + bias2
return (scvmat)
}
###############################################################################
## Find the bandwidth that minimises the SCV for 1 to 6 dim
##
## Parameters
## x - data values
## pre - "scale" - pre-scaled data
## - "sphere"- pre-sphered data
## Hstart - initial bandwidth matrix
##
## Returns
## H_SCV
###############################################################################
hscv <- function(x, nstage=2, binned=TRUE, bgridsize, plot=FALSE)
{
sigma <- sd(x)
n <- length(x)
d <- 1
hnorm <- sqrt((4/(n*(d + 2)))^(2/(d + 4)) * var(x))
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
hmin <- 0.1*hnorm
hmax <- 2*hnorm
bin.par <- binning(x=x, bgridsize=bgridsize, h=hnorm)
if (nstage==1)
{
psihat6 <- psins.1d(r=6, sigma=sigma)
psihat10 <- psins.1d(r=10, sigma=sigma)
}
else if (nstage==2)
{
g1 <- (2/(7*n))^(1/9)*2^(1/2)*sigma
g2 <- (2/(11*n))^(1/13)*2^(1/2)*sigma
psihat6 <- kfe.1d(x=x, bin.par=bin.par, binned=binned, deriv.order=6, g=g1, inc=1)
psihat10 <- kfe.1d(x=x, bin.par=bin.par, binned=binned, deriv.order=10, g=g2, inc=1)
}
g3 <- (-6/((2*pi)^(1/2)*psihat6*n))^(1/7)
g4 <- (-210/((2*pi)^(1/2)*psihat10*n))^(1/11)
psihat4 <- kfe.1d(x=x, bin.par=bin.par, binned=binned, deriv.order=4, g=g3, inc=1)
psihat8 <- kfe.1d(x=x, bin.par=bin.par, binned=binned, deriv.order=8, g=g4, inc=1)
C <- (441/(64*pi))^(1/18) * (4*pi)^(-1/5) * psihat4^(-2/5) * psihat8^(-1/9)
scv.1d.temp <- function(h)
{
return(scv.1d(x=x, bin.par=bin.par, h=h, g=C*n^(-23/45)*h^(-2), binned=binned, inc=1))
}
if (plot)
{
hseq <- seq(hmin,hmax, length=400)
hscv.seq <- rep(0, length=length(hseq))
for (i in 1:length(hseq))
hscv.seq[i] <- scv.1d.temp(hseq[i])
plot(hseq, hscv.seq, type="l", xlab="h", ylab="SCV(h)")
}
opt <- optimise(f=scv.1d.temp, interval=c(hmin, hmax))$minimum
if (n >= 1e5) warning("hscv is not always stable for large samples")
return(opt)
}
Hscv <- function(x, nstage=2, pre="sphere", pilot, Hstart, binned, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim")
{
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d > 4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
if(!is.matrix(x)) x <- as.matrix(x)
if (missing(pilot)) { if (d==2 & r==0) pilot <- "samse" else pilot <- "dscalar" }
pilot1 <- match.arg(pilot, c("amse", "samse", "unconstr", "dunconstr", "dscalar"))
pre1 <- match.arg(pre, c("scale", "sphere"))
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
if (pilot1=="amse" & (d>2 | r>0)) stop("amse pilot selectors not defined for d>2 and/or r>0")
if ((pilot1=="samse" | pilot1=="unconstr") & r>0) stop("dscalar or dunconstr pilot selectors are better for deriv.order>0")
if (pre1=="scale")
{
x.star <- pre.scale(x)
S12 <- diag(sqrt(diag(var(x))))
Sinv12 <- chol2inv(chol(S12))
}
else if (pre1=="sphere")
{
x.star <- pre.sphere(x)
S12 <- matrix.sqrt(var(x))
Sinv12 <- chol2inv(chol(S12))
}
RK <- (4*pi)^(-d/2)
if (binned)
{
if (pilot1=="unconstr" | pilot1=="dunconstr")
H.max <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x)
else
H.max <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
if (default.bflag(d=d, n=n))
bin.par <- binning(x=x.star, bgridsize=bgridsize, H=matrix.sqrt(H.max))
}
if (pilot1=="unconstr")
{
## Gu pilot matrix is on data scale
Gu <- Gunconstr.scv(x=x, binned=binned, bgridsize=bgridsize, verbose=verbose, nstage=nstage-1, optim.fun=optim.fun)
if (missing(Hstart)) Hstart <- Hns(x=x, deriv.order=r)
}
else if (pilot1=="dunconstr")
{
## Gu pilot matrix is on data scale
Gu <- Gdunconstr(x=x, d=d, r=r, n=n, nstage=nstage, verbose=verbose, binned=binned, scv=TRUE, optim.fun=optim.fun)
if (missing(Hstart)) Hstart <- Hns(x=x, deriv.order=r)
}
else if (pilot1=="dscalar")
{
## Gs is on pre-transformed data scale
g2r4 <- gdscalar(x=x.star, d=d, r=r, n=n, nstage=nstage, verbose=verbose, scv=TRUE, binned=binned)
Gs <- g2r4^2*diag(d)
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
}
else
{
## Gs is on transformed data scale
Hamise <- Hpi(x=x.star, nstage=1, deriv.order=r, pilot=pilot, pre="sphere", binned=TRUE, bgridsize=bgridsize, verbose=verbose, optim.fun=optim.fun)
if (any(is.na(Hamise)))
{
warning("Pilot bandwidth matrix is NA - replaced with maximally smoothed")
Hamise <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
}
gs <- gamse.scv(x.star=x.star, d=d, Sigma.star=var(x.star), Hamise=Hamise, n=n, binned=binned, bgridsize=bgridsize, verbose=verbose, nstage=nstage-1)
Gs <- gs^2*diag(d)
## use normal reference bandwidth as initial condition
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
else Hstart <- Sinv12 %*% Hstart %*% Sinv12
}
## SCV is estimate of AMISE
scv.mat.temp <- function(vechH)
{
H <- invvech(vechH) %*% invvech(vechH)
if (pilot1=="samse" | pilot1=="amse" | pilot1=="dscalar") { Gpilot <- Gs; xx <- x.star }
else if (pilot1=="unconstr" | pilot1=="dunconstr") { Gpilot <- Gu; xx <- x }
if (default.bflag(d=d, n=n))
scvm <- scv.mat(x=xx, H=H, G=Gpilot, binned=binned, bin.par=bin.par, verbose=FALSE, deriv.order=r)
else
scvm <- scv.mat(x=xx, H=H, G=Gpilot, binned=binned, verbose=FALSE, deriv.order=r)
return(scvm)
}
Hstart <- matrix.sqrt(Hstart)
if (optim.fun=="nlm")
{
result <- nlm(p=vech(Hstart), f=scv.mat.temp, print.level=2*as.numeric(verbose))
H <- invvech(result$estimate) %*% invvech(result$estimate)
amise.star <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(vech(Hstart), scv.mat.temp, method="BFGS", control=list(trace=as.numeric(verbose), REPORT=1))
H <- invvech(result$par) %*% invvech(result$par)
amise.star <- result$value
}
if (!(pilot1 %in% c("dunconstr","unconstr"))) H <- S12 %*% H %*% S12 ## back-transform
if (!amise) return(H)
else return(list(H = H, SCV.star=amise.star))
}
Hscv.diag <- function(x, nstage=2, pre="scale", pilot, Hstart, binned, bgridsize, amise=FALSE, deriv.order=0, verbose=FALSE, optim.fun="optim")
{
n <- nrow(x)
d <- ncol(x)
r <- deriv.order
RK <- (4*pi)^(-d/2)
if (missing(binned)) binned <- default.bflag(d=d,n=n)
if (d > 4) binned <- FALSE
if (missing(bgridsize)) bgridsize <- default.bgridsize(d)
if(!is.matrix(x)) x <- as.matrix(x)
if (missing(pilot)) { if (d==2 & r==0) pilot <- "samse" else pilot <- "dscalar" }
pilot1 <- match.arg(pilot, c("amse", "samse", "unconstr", "dunconstr", "dscalar"))
pre1 <- match.arg(pre, c("scale", "sphere"))
optim.fun <- match.arg(optim.fun, c("nlm", "optim"))
if (pilot1=="amse" & (d>2 | r>0)) stop("samse pilot selectors are better for higher dimensions and/or deriv.order>0")
if (pilot1=="samse" & r>0) stop("dscalar pilot selectors are better for deriv.order>0")
if (pilot1=="unconstr" | pilot1=="dunconstr") stop("Unconstrained pilot selectors are not suitable for Hscv.diag")
if (pre1=="sphere") stop("Using pre-sphering doesn't give a diagonal bandwidth matrix")
if (pre1=="scale")
{
x.star <- pre.scale(x)
S12 <- diag(sqrt(diag(var(x))))
Sinv12 <- chol2inv(chol(S12))
}
else if (pre1=="sphere")
{
x.star <- pre.sphere(x)
S12 <- matrix.sqrt(var(x))
Sinv12 <- chol2inv(chol(S12))
}
if (binned)
{
H.max <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
bin.par.star <- binning(x=x.star, bgridsize=bgridsize, H=H.max)
}
if (pilot1=="dscalar")
{
## Gs is on pre-transformed data scale
g2r4 <- gdscalar(x=x.star, r=r, n=n, d=d, verbose=verbose, nstage=nstage, scv=TRUE, binned=binned)
Gs <- g2r4^2*diag(d)
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
}
else
{
## Gs is on transformed data scale
Hamise <- Hpi(x=x.star, nstage=1, pilot=pilot, pre="sphere", binned=binned, bgridsize=bgridsize, verbose=verbose, optim.fun=optim.fun)
if (any(is.na(Hamise)))
{
warning("Pilot bandwidth matrix is NA - replaced with maximally smoothed")
Hamise <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)
}
gs <- gamse.scv(x.star=x.star, d=d, Sigma.star=var(x.star), Hamise=Hamise, n=n, binned=binned, bgridsize=bgridsize, verbose=verbose, nstage=nstage-1)
Gs <- gs^2*diag(d)
## use normal reference bandwidth as initial condition
if (missing(Hstart)) Hstart <- Hns(x=x.star, deriv.order=r)
else Hstart <- Sinv12 %*% Hstart %*% Sinv12
}
scv.mat.temp <- function(diagH)
{
H <- diag(diagH) %*% diag(diagH)
if (default.bflag(d=d, n=n)) scvm <- scv.mat(x.star, H, Gs, binned=binned, verbose=FALSE, bin.par=bin.par.star, deriv.order=r)
else scvm <- scv.mat(x.star, H, Gs, binned=binned, verbose=FALSE, deriv.order=r)
return(scvm)
}
Hstart <- matrix.sqrt(Hstart)
if (optim.fun=="nlm")
{
result <- nlm(p=diag(Hstart), f=scv.mat.temp, print.level=2*as.numeric(verbose))
H <- diag(result$estimate) %*% diag(result$estimate)
amise.star <- result$minimum
}
else if (optim.fun=="optim")
{
result <- optim(diag(Hstart), scv.mat.temp, method="Nelder-Mead", control=list(trace=as.numeric(verbose), REPORT=1))
H <- diag(result$par) %*% diag(result$par)
amise.star <- result$value
}
## back-transform
H <- S12 %*% H %*% S12
if (!amise) return(H)
else return(list(H = H, SCV.star=amise.star))
}
##############################################################################
## Normal scale selector H_ns for kernel density derivate estimators
##############################################################################
Hns <- function(x, deriv.order=0)
{
if (is.vector(x)) { n <- 1; d <- length(x) }
else { n <- nrow(x); d <- ncol(x) }
r <- deriv.order
H <- (4/(n*(d+2*r+2)))^(2/(d+2*r+4)) * var(x)
return(H)
}
Hns.diag <- function(x)
{
if (is.vector(x)) { n <- 1; d <- length(x) }
else { n <- nrow(x); d <- ncol(x) }
S <- var(x)
Delta <- chol2inv(chol(diag(diag(S))))%*%S
svD <- svd(Delta)
Deltainv <- svD$v %*% diag(1/svD$d) %*% t(svD$u)
H <- ((4*d*det(Delta)^(1/2))/(2*tr(Deltainv%*%Deltainv) + (tr(Deltainv))^2))^(2/(d+4))*diag(diag(S))*n^(-2/(d+4))
return(H)
}
hns <- function(x, deriv.order=0)
{
n <- length(x)
d <- 1
r <- deriv.order
h <- (4/(n*(d+2*r+2)))^(1/(d+2*r+4))*sd(x)
return(h)
}
#######################################################################
## Normal scale G_ns for kernel functional estimators
## scv = flag for SCV bias annihiliation constant
#######################################################################
Gns <- function(r,n,Sigma,scv=FALSE)
{
d <- ncol(Sigma)
const <- ifelse(scv,1,2)
G <- (2/((n*(d+r))))^(2/(d+r+2))*const*Sigma
return(G)
}
Gns.search <- function(G, f, n=10)
{
Gstart.vec <- numeric(n)
for (i in 1:length(Gstart.vec))
{
Gstart <- matrix.sqrt(i/length(Gstart.vec)*G)
Gstart.vec[i] <- f(vech(Gstart))
}
i <- which.min(Gstart.vec)/length(Gstart.vec)
Gstart <- i*G
return(Gstart)
}
##############################################################################
## Normal mixture selector
##############################################################################
Hnm <- function(x, deriv.order=0, G=1:9, subset.ind, mise.flag=FALSE, verbose=FALSE, ...)
{
if (!requireNamespace("mclust", quietly=TRUE)) stop("Install the mclust package as it is required.", call.=FALSE)
modelNames <- "VVV"
if (!missing(subset.ind)) nmixt.fit <- mclust::Mclust(x[subset.ind,], G=G, verbose=verbose, modelNames=modelNames, ...)
else nmixt.fit <- mclust::Mclust(x, G=G, verbose=verbose, modelNames=modelNames, ...)
if (is.vector(x)) { d <- length(x); n <- 1 } else { d <- ncol(x); n <- nrow(x) }
mus <- t(nmixt.fit$parameters$mean)
Sigmas <- matrix(nmixt.fit$parameters$variance$sigma, byrow=TRUE, ncol=d)
props <- nmixt.fit$parameters$pro
if (mise.flag) H.nm <- Hmise.mixt(samp=n, mus=mus, Sigmas=Sigmas, props=props, deriv.order=deriv.order)
else H.nm <- Hamise.mixt(samp=n, mus=mus, Sigmas=Sigmas, props=props, deriv.order=deriv.order)
return(H.nm)
}
Hnm.diag <- function(x, deriv.order=0, G=1:9, subset.ind, mise.flag=FALSE, verbose=FALSE, ...)
{
if (!requireNamespace("mclust", quietly=TRUE)) stop("Install the mclust package as it is required.", call.=FALSE)
modelNames <- "VVI"
if (!missing(subset.ind)) nmixt.fit <- mclust::Mclust(x[subset.ind,], G=G, verbose=verbose, modelNames=modelNames, ...)
else nmixt.fit <- mclust::Mclust(x, G=G, verbose=verbose, modelNames=modelNames, ...)
if (is.vector(x)) { d <- length(x); n <- 1 } else { d <- ncol(x); n <- nrow(x) }
mus <- t(nmixt.fit$parameters$mean)
Sigmas <- matrix(nmixt.fit$parameters$variance$sigma, byrow=TRUE, ncol=d)
props <- nmixt.fit$parameters$pro
if (mise.flag) H.nm <- Hmise.mixt.diag(samp=n, mus=mus, Sigmas=Sigmas, props=props, deriv.order=deriv.order)
else H.nm <- Hamise.mixt.diag(samp=n, mus=mus, Sigmas=Sigmas, props=props, deriv.order=deriv.order)
return(H.nm)
}
hnm <- function(x, deriv.order=0, G=1:9, subset.ind, mise.flag=FALSE, verbose=FALSE, ...)
{
if (!missing(subset.ind)) nmixt.fit <- mclust::Mclust(x[subset.ind], G=G, verbose=verbose, ...)
else nmixt.fit <- mclust::Mclust(x, G=G, verbose=verbose, ...)
mus <- nmixt.fit$parameters$mean
sigmas <- sqrt(nmixt.fit$parameters$variance$sigma)
props <- nmixt.fit$parameters$pro
n <- length(x)
if (mise.flag) h.nm <- hmise.mixt(samp=n, mus=mus, sigmas=sigmas, props=props, deriv.order=deriv.order)
else h.nm <- hamise.mixt(samp=n, mus=mus, sigmas=sigmas, props=props, deriv.order=deriv.order)
return(h.nm)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.