Nothing
# kernels.R
# - useful kernels...
# NOTES:
# - list should increase..
# - 'gaussK' is very slowly to do computations, better a compact supported kernel...
# ERRORS:
#
# Kernels (ok)
#
# Gaussian kernel
gaussK <- function(x) dnorm(x,0,1)
## OK
attr(gaussK,"RK") <- 1/( 2 * sqrt(pi) )
attr(gaussK,"RdK") <- 1/( 4 * sqrt(pi) )
attr(gaussK,"mu0K") <- 1
attr(gaussK,"mu2K") <- 1
attr(gaussK,"K4") <- 0.1994711
attr(gaussK,"dom") <- c(-Inf,Inf)
# Gaussian kernel(locfit version)
gaussKlf <- function(x) exp(-(2.25 * x)^2);
## OK
attr(gaussKlf,"RK") <- 0.557029
attr(gaussKlf,"RdK") <- 2.81996
attr(gaussKlf,"mu0K") <- 0.787757
attr(gaussKlf,"mu2K") <- 0.0778032
attr(gaussKlf,"K4") <- 0.2444259
attr(gaussKlf,"dom") <- c(-Inf,Inf)
# Squared kernel
SqK <- function(x) ifelse( abs(x) <= 1 , 1 , 0 );
## OK
attr(SqK,"RK") <- 2.0
attr(SqK,"RdK") <- 0.0
attr(SqK,"mu0K") <- 2.0
attr(SqK,"mu2K") <- 0.666667
attr(SqK,"K4") <- NA
# Triangle kernel
TrianK <- function(x) ifelse( abs(x) <= 1 , (1 - abs(x)) , 0 );
## OK
attr(TrianK,"RK") <- 0.666667
attr(TrianK,"RdK") <- 2
attr(TrianK,"mu0K") <- 1.
attr(TrianK,"mu2K") <- 0.166667
attr(TrianK,"K4") <- 0.4793729
attr(TrianK,"dom") <- c(-1,1)
# Epanechnikov kernel
EpaK <- function(x) ifelse( abs(x) <= 1 , 3/4*(1 - x^2) , 0 );
## OK
attr(EpaK,"RK") <- 0.6
attr(EpaK,"RdK") <- 1.5
attr(EpaK,"mu0K") <- 1.
attr(EpaK,"mu2K") <- 0.2
attr(EpaK,"K4") <- 0.4337657
attr(EpaK,"dom") <- c(-1,1)
# Epanechnikov kernel
Epa2K <- function(x) ifelse( abs(x) <= 1 , 3/(4*sqrt(5))*(1 - x^2/5) , 0 );
## OK
attr(Epa2K,"RK") <- 0.1968
attr(Epa2K,"RdK") <- 0.012
attr(Epa2K,"mu0K") <- 0.626099
attr(Epa2K,"mu2K") <- 0.196774
attr(Epa2K,"K4") <- NA
attr(Epa2K,"dom") <- c(-1,1)
# Quartic kernel
QuartK <- function(x) ifelse( abs(x) <= 1 , 15/16*(1 - x^2)^2 , 0 );
## OK
attr(QuartK,"RK") <- 0.714286
attr(QuartK,"RdK") <- 2.14286
attr(QuartK,"mu0K") <- 1.
attr(QuartK,"mu2K") <- 0.142857
attr(QuartK,"K4") <- 0.5164128
attr(QuartK,"dom") <- c(-1,1)
biweigK <- function(x) ifelse( abs(x) <= 1 , (1 - x^2)^2 , 0 );
## OK
attr(biweigK,"RK") <- 0.812698;
attr(biweigK,"RdK") <- 2.4381;
attr(biweigK,"mu0K") <- 1.06667;
attr(biweigK,"mu2K") <- 0.152381;
attr(biweigK,"K4") <- 0.6685162
attr(biweigK,"dom") <- c(-1,1)
# Triweigth kernel
TriweigK <- function(x) ifelse( abs(x) <= 1 , 35/32*(1 - x^2)^3 , 0 );
## OK
attr(TriweigK,"RK") <- 0.815851
attr(TriweigK,"RdK") <- 3.18182
attr(TriweigK,"mu0K") <- 1.
attr(TriweigK,"mu2K") <- 0.111111
attr(TriweigK,"K4") <- 0.5879012
attr(TriweigK,"dom") <- c(-1,1)
tricubK <- function(x) ifelse( abs(x) <= 1 , 70/81*(1 - abs(x)^3)^3 , 0 );
## OK
attr(tricubK,"RK") <- 0.708502;
attr(tricubK,"RdK") <- 2.24599;
attr(tricubK,"mu0K") <- 1;
attr(tricubK,"mu2K") <- 0.144033;
attr(tricubK,"K4") <- 0.5879012
attr(tricubK,"dom") <- c(-1,1)
tricubKlf <- function(x) ifelse( abs(x) <= 1 , (1 - abs(x)^3)^3 , 0 );
## OK, (locfit version)
attr(tricubKlf,"RK") <- 0.94867;
attr(tricubKlf,"RdK") <- 3.00733
attr(tricubKlf,"mu0K") <- 1.15714;
attr(tricubKlf,"mu2K") <- 0.166667;
attr(tricubKlf,"K4") <- NA
attr(tricubKlf,"dom") <- c(-1,1)
# cosin kernel
CosK <- function(x) ifelse( abs(x) <= 1 , pi/4*cos(pi*x/2) , 0 );
## OK
attr(CosK,"RK") <- 0.6168
attr(CosK,"RdK") <- 1.52202
attr(CosK,"mu0K") <- 1.
attr(CosK,"mu2K") <- 0.189431
attr(CosK,"K4") <- 0.4464387
attr(CosK,"dom") <- c(-1,1)
.kernelList <- c( "gaussK", "EpaK", "Epa2K", "TrianK", "QuartK", "biweigK",
"TriweigK", "tricubK","CosK", "SqK" )
# R(K) wrapper
RK <- function(K) return( attr(K,"RK") )
# R(K') wrapper
RdK <- function(K) return( attr(K,"RdK") )
# mu2K wrapper
mu2K <- function(K) return( attr(K,"mu2K") )
# mu2K wrapper
mu0K <- function(K) return( attr(K,"mu0K") )
# K4 wrapper
K4 <- function(K) return( attr(K,"K4") )
# dom wrapper
dom <- function(K) return( attr(K,"dom") )
##
## Building equivalent kernel
##
equivKernel <- function(kernel,nu,deg,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##OK
# p = local polinomial order.
# nu = which of the equivalent kernels do you want.
# kernel = kernel.
# NOTES:
# Al menos para p=1,1 funciona
{
mu <- array( dim=c(deg+1,deg+1) )
for(i in 0:deg)
for(j in 0:deg)
mu[i+1,j+1] <- computeMu(i+j,kernel,lower,upper,subdivisions)
invMu <- solve(mu)
res <- function(t)
{
return( (invMu %*% outer(0:deg,t,function(x,y) y^x))[nu+1,]*kernel(t) )
##return( (invMu %*% outer(nu+1,t,function(x,y) y^x))*kernel(t) )
}
return( res )
}
##
## Building C_{\nu p}(K)
##
cteNuK <- function(nu,p,kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
# nu = ;
# p = ;
# kernel = ;
# NOTES:
# Al menos para EpaK, y gaussK y TriweigK con p=0,1,... la cosa funciona...
# Es decir funciona...
{
if( p==0 || p==1 ){
a <- RK(kernel)
b <- mu2K(kernel)
return( ((p+1)^2*(2*nu+1)*a/(2*(p+1-nu)*b^2))^(1/(2*p+3)) )
}else{
f <- 1
for(i in 2:(p+1))
f <- f * i
eqKern <- equivKernel(kernel,nu,p,lower,upper,subdivisions)
a <- computeRK(eqKern,lower,upper,subdivisions)
b <- computeMu(p+1,eqKern,lower,upper,subdivisions)
return( (f^2*(2*nu+1)*a/(2*(p+1-nu)*b^2))^(1/(2*p+3)) )
}
}
##
## Building adj_{\nu p}(K), to compute derivative estimation bandwidth.
##
adjNuK <- function(nu,p,kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
# nu = ;
# p = ;
# kernel = ;
# RETURNS
# Just 'cteNuK(nu,p,kernel)/cteNuK(0,p,kernel)', useful for derivative
# bandwidth computations.
{
a <- cteNuK(nu,p,kernel,lower,upper,subdivisions)
b <- cteNuK(0,p,kernel,lower,upper,subdivisions)
return( a/b )
}
##
## Some kernel miscelanea functions
##
Kconvol <- function(kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
# RETURNS
# given Kernel Convolution function. Computed by Num. Int.
# kernel = Kernel use in computation
# order = currently onder two is the only one computed...
# NOTES:
# Use 'integrate' from R, seems to be quick and accurate!!!
{
ff <- function(x) integrate(
function(u) return( kernel(u) * kernel(x-u) ),
lower,upper,subdivisions=subdivisions)$value
res <- function(x) sapply(x,ff)
return( res )
}
computeRK <- function(kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
#
# kernel = Kernel use in computation
# lower = lower limit of interval.
# upper = upper limit of interval.
# subdivisions = number of subdivisions used to compute integral.
# Computes RK for the given kernel.
{
f <- function(u) kernel(u)^2
return( integrate(f,lower,upper,subdivisions=subdivisions)$value )
}
computeK4 <- function(kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
#
# kernel = Kernel use in computation
# lower = lower limit of interval.
# upper = upper limit of interval.
# subdivisions = number of subdivisions used to compute integral.
# Computes K4 for the given kernel, the autoconvolution fo the
# autoconvolution for the given kernel .
{
f <- Kconvol(kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
ff <- function(x) f(x)^2
return( integrate(ff,2*lower,2*upper,subdivisions=subdivisions)$value )
}
# computeRdK <- function(kernel,lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
# {
# dd <-
# f <- function(u) kernel(u)^2
# return( integrate(f,lower,upper,subdivisions)$value )
# }
computeMu0 <- function(kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
#
# kernel = Kernel use in computation
# lower = lower limit of interval.
# upper = upper limit of interval.
# subdivisions = number of subdivisions used to compute integral.
# Computes mu_0 for the given kernel.
{
return( integrate(kernel,lower,upper,subdivisions=subdivisions)$value )
}
computeMu <- function(i,kernel,
lower=dom(kernel)[[1]],upper=dom(kernel)[[2]],subdivisions=25)
##
#
# kernel = Kernel use in computation
# lower = lower limit of interval.
# upper = upper limit of interval.
# subdivisions = number of subdivisions used to compute integral.
# Computes mu_i for the given kernel.
{
f <- function(u) u^i*kernel(u)
return( integrate(f,lower,upper,subdivisions=subdivisions)$value )
}
# locPolWeights.R
#
# NOTAS:
# ERRORES:
locPolWeights <- function (x, xeval, deg, bw, kernel,
weig=rep(1,length(x)))
# Suggestion from munevvere@hotmail.com
{
stopifnot( length(x)==length(weig), bw>0, deg>=0 )
den <- array( dim=length(xeval),dimnames=list(x=1:length(xeval)) )
res <- array( dim=c(length(xeval),deg+1,length(x)),
dimnames=list(x=1:length(xeval),deg=0:deg,xData=1:length(x)))
for(i in 1:length(xeval))
{
xx <- xeval[i]
dMat <- (x-xx)/bw
xx <- outer(dMat, 0:deg, function(a, b) a^b)
w <- diag(kernel(dMat)*weig)
aux <- t(xx) %*% w
sMat <- aux %*% xx
den[i] <- det(sMat)
res[i,,] <- solve(sMat,aux)/bw^(0:deg)
}
invisible( list(den=den*bw^deg,locWeig=res[,1,],allWeig=res) )
}
#
# locpol.R
# Interface para el polinomio local en C
# NOTAS:
# - A?adir todos los kernels...
# - Ojito, devuelve 0 si det(X^TWX)=0!!
# ERRORES:
## old s3 classes stuff
# .First.lib <- function(lib, pkg) {
# library.dynam("locpol", pkg, lib)
# }
## s4 classes dim lib load stuff
.onLoad <- function(lib, pkg) {
library.dynam("locpol", pkg, lib)
}
.maxEvalPts <- 5000
.lokestOptInt <- c(0.0005,1.5)
selKernel <- function(kernel)
{
if( RK(kernel)==0.6 ) ## Epanechnikov kernel(ok R)
return(1)
else if ( RK(kernel)==0.1968 ) ## 2nd. Epanechnikov kernel(ok R)
return(2)
else if ( RK(kernel)==0.666667 ) ## Triangle kernel(ok R)
return(3)
else if ( RK(kernel)==0.714286 ) ## Quartic kernel(ok R)
return(4)
else if ( RK(kernel)==0.812698 ) ## biweight kernel(ok R)
return(5)
else if ( RK(kernel)==0.815851 ) ## Triweigth kernel(ok R)
return(6)
else if ( RK(kernel)==0.708502 ) ## tricube kernel(ok R)
return(7)
else if ( RK(kernel)==0.6168 ) ## cosin kernel(ok R)
return(9)
else if ( RK(kernel)==2.0 ) ## Square kernel(ok R)
return(10)
else ## gaussian kernel
return(0)
}
locWeightsEval <- function(lpweig,y)
##OK
#
# lpweig = an x-data frame.
# y = eval. points.
{
stopifnot( ncol(lpweig) == length(y), nrow(lpweig)<.maxEvalPts )
return( lpweig %*% y )
}
locWeightsEvalC <- function(lpweig,y)
##OK
#
# lpweig = an x-data frame.
# y = eval. points.
{
stopifnot( ncol(lpweig) == length(y), nrow(lpweig)<.maxEvalPts )
res <- .C("locWeightsEvalxx",
as.double(lpweig),as.integer(nrow(lpweig)),
as.double(y),as.integer(length(y)),
res=double(nrow(lpweig)), PACKAGE="locpol"
)
#rownames(res$res) <- rownames(lpweig)
return( res$res )
}
simpleSmootherC <- function(x,y,xeval,bw,kernel,weig = rep(1,length(y)))
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'simpleSmoother'.
# Computes m(x)f(x)mu_0.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0)
Ktype <- selKernel(kernel)
reg <- .C("simpleSmoother",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(y),as.double(weig),
as.integer(length(y)),as.double(bw), as.integer(Ktype),
res = double(length(xeval)), PACKAGE="locpol"
)$res
res <- data.frame(x = xeval, reg = reg)
return( res )
}
PRDenEstC <- function(x,xeval,bw,kernel,weig = rep(1,length(x)))
##OK
#
# x = an x vector.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'parzenRossen'.
# Computes the density of x
{
stopifnot(length(x)==length(weig), bw > 0)
Ktype <- selKernel(kernel)
den <- .C("parzenRossen",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(weig),as.integer(length(x)),
as.double(bw), as.integer(Ktype),
res = double(length(xeval)), PACKAGE="locpol"
)$res
res <- data.frame(x = xeval, den = den)
return( res )
}
locCteSmootherC <- function(x,y,xeval,bw,kernel,weig = rep(1,length(y)))
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'locCteSmoother'.
# Computes m(x), Nadaraya-Watson estimator, as a by--product we
# got the marginal density for X.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0,
length(xeval)<.maxEvalPts)
Ktype <- selKernel(kernel)
res <- .C("locCteSmoother",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(y),as.double(weig),
as.integer(length(y)),as.double(bw), as.integer(Ktype),
den=double(length(xeval)),res=double(length(xeval)), PACKAGE="locpol"
)
res <- data.frame(x = xeval, beta0 = res$res, den=res$den)
return( res )
}
locCteWeightsC <- function(x,xeval,bw,kernel,weig = rep(1,length(x)))
##OK
#
# x = x vector.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'locCteSmoother'.
# Computes m(x), Nadaraya-Watson estimator, as a by--product we
# got the marginal density for X.
{
stopifnot(length(x)==length(weig), bw > 0, length(xeval)<.maxEvalPts)
Ktype <- selKernel(kernel)
res <- .C("locCteWeights",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(weig),as.integer(length(x)),
as.double(bw), as.integer(Ktype),
den=double(length(xeval)),res=mat.or.vec(length(xeval),length(x)),
PACKAGE="locpol"
)
rownames(res$res) <- xeval
colnames(res$res) <- 1:length(x)
return( list(den=res$den,locWeig=res$res) )
}
locLinSmootherC <- function(x,y,xeval,bw,kernel,weig = rep(1,length(y)))
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'locLinSmoother'.
# Computes m(x), Local linear estimator, as a by--product we
# got f(x)^2, the squared marginal density for X.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0,
length(xeval)<.maxEvalPts)
Ktype <- selKernel(kernel)
res <- .C("locLinSmoother",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(y),as.double(weig),
as.integer(length(y)),
as.double(bw), as.integer(Ktype),den=double(length(xeval)),
beta0=double(length(xeval)),beta1=double(length(xeval)),
PACKAGE="locpol"
)
res <- data.frame(x=xeval, beta0=res$beta0, beta1=res$beta1, den=res$den)
return( res )
}
locLinWeightsC <- function(x,xeval,bw,kernel,weig = rep(1,length(x)))
##OK
#
# x = x vector.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'locCteSmoother'.
# Computes m(x), Nadaraya-Watson estimator, as a by--product we
# got the marginal density for X.
{
stopifnot(length(x)==length(weig), bw > 0, length(xeval)<.maxEvalPts)
Ktype <- selKernel(kernel)
res <- .C("locLinWeights",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(weig),
as.integer(length(x)),
as.double(bw), as.integer(Ktype),
den=double(length(xeval)),res=mat.or.vec(length(xeval),length(x)),
PACKAGE="locpol"
)
rownames(res$res) <- xeval
colnames(res$res) <- 1:length(x)
return( list(den=res$den,locWeig=res$res) )
}
locCuadSmootherC <- function(x,y,xeval,bw,kernel,weig = rep(1,length(y)))
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# An interface to C function 'locCuadSmoother'.
# Computes m(x), Local cuadratic estimator, as a by--product we
# got f(x)^3, the squared marginal density for X.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0,
length(xeval)<.maxEvalPts)
Ktype <- selKernel(kernel)
res <- .C("locCuadSmoother",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(y),as.double(weig),
as.integer(length(y)),as.double(bw), as.integer(Ktype),
den=double(length(xeval)), beta0=double(length(xeval)),
beta1=double(length(xeval)), beta2=double(length(xeval)),
PACKAGE="locpol"
)
res <- data.frame(x = xeval, beta0=res$beta0, beta1=res$beta1,
beta2=res$beta2, den=res$den)
return( res )
}
locPolSmootherC <- function(x,y,xeval,bw,deg,kernel,DET=FALSE,
weig=rep(1,length(y)))
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# DET = flag to compute the determinant of t(X)WX(X local desing matrix)
# An interface to C function 'locPolSmootherC'.
# Computes m(x), Local pol. estimator, as a by--product we
# got f(x)^(p+1), the squared marginal density for X.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0,
length(xeval)< .maxEvalPts,0<=deg && deg<10)
Ktype <- selKernel(kernel)
res <- .C("locPolSmoother",
as.double(xeval), as.integer(length(xeval)),
as.double(x), as.double(y), as.double(weig),
as.integer(length(y)), as.double(bw), as.integer(deg),
as.integer(Ktype), as.integer(DET),
den=double(length(xeval)),
beta=mat.or.vec(length(xeval),(deg+1)), PACKAGE="locpol"
)
res <- data.frame(x=xeval, beta=matrix(res$beta,ncol=deg+1),res$den)
colnames(res) <- c("x",paste("beta",0:deg,sep=""),"den")
if( !DET ) res$den <- NULL
return( res )
}
simpleSqSmootherC <- function(x,y,xeval,bw,kernel)
##OK
#
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# kernel = kernel used in estimation.
# An interface to C function 'simpleSmoother'. It returns
# E[y^2|x]f(x)R(K), can be useful to estimate variance...
{
stopifnot(length(x)==length(y), bw > 0)
Ktype <- selKernel(kernel)
reg <- .C("simpleSqSmoother",
as.double(xeval),as.integer(length(xeval)),
as.double(x),as.double(y),
as.integer(length(y)),
as.double(bw), as.integer(Ktype),
res = double(length(xeval)), PACKAGE="locpol")$res
res <- data.frame(x = xeval, reg = reg)
return( res )
}
##
## Leave One out computations.
##
looLocPolSmootherC <- function(x,y,bw,deg,kernel,weig=rep(1,length(y)),DET=FALSE)
##OK
# x, y = x and y vectors.
# xeval = eval. points.
# bw = bandwidth.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# DET = flag to compute the determinant of t(X)WX(X local desing matrix)
# An interface to C function 'looLocPolSmootherC'. Estimation of
# m(x_i) i=1,...n without ith observation (x_i,y_i).
# Computes m(x), Local pol. estimator, as a by--product we
# got f(x)^(p+1), the squared marginal density for X.
{
stopifnot(length(x)==length(y), length(y)==length(weig), bw > 0,
length(x)<.maxEvalPts,0<=deg && deg<10)
Ktype <- selKernel(kernel)
res <- .C("looLocPolSmoother",
as.double(x), as.double(y), as.double(weig),
as.integer(length(y)), as.double(bw), as.integer(deg),
as.integer(Ktype), as.integer(DET),
den=double(length(y)), beta=mat.or.vec(length(y),(deg+1)),
PACKAGE="locpol"
)
res <- data.frame(x=x, beta=matrix(res$beta,ncol=deg+1),res$den)
colnames(res) <- c("x",paste("beta",0:deg,sep=""),"den")
return( res )
}
##
## Bandwidth selection procedures...
##
denCVBwSelC <- function(x,kernel=gaussK,weig=rep(1,length(x)),
interval=.lokestOptInt)
##OK
#
# x = x vector
# xeval = eval. points.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# Computes Cross Validation bandwidth selector for the
# density estimator...
{
stopifnot(length(x)==length(weig))
Ktype <- selKernel(kernel)
cvFunc <- function(h) .C("denCVBwEval",as.double(h),
as.double(x),as.double(weig),as.integer(length(x)),
as.integer(Ktype),res = double(1), PACKAGE="locpol")$res
return( optimise(cvFunc,interval)$minimum )
}
regCVBwSelC <- function(x,y,deg,kernel=gaussK,weig=rep(1,length(y)),
interval=.lokestOptInt)
##OK
#
# dat = an x-data frame.
# xeval = eval. points.
# weig = vector of weigths for observations.
# kernel = kernel used in estimation.
# Computes Cross Validation bandwidth selector for the
# regression function estimator...
{
stopifnot(length(x)==length(y), length(y)==length(weig), 0<=deg && deg<10 )
Ktype <- selKernel(kernel)
cvFunc <- function(h) .C("regCVBwEvalB",as.double(h),
as.double(x), as.double(y), as.double(weig),
as.integer(length(y)), as.integer(deg),
as.integer(Ktype),res = double(1), PACKAGE="locpol")$res
return( optimise(cvFunc,interval)$minimum )
}
# locpolClass.R
# local polynomial class
# NOTES:
# -- Bias estimations !! ??
# -- Keep track of errors... So you'll be able to give always
# a suggestion in case errors appears
# -- Imporve error messages !!
# ERRORS:
locpol <- function(formula, data, weig=rep(1,nrow(data)),
bw=NULL, kernel=EpaK, deg=1, xeval=NULL,xevalLen=100)
##
{
## checking
stopifnot(nrow(data)==length(weig))
## compute result
res <- list()
res$bw <- bw
res$KName <- match.call()
res$kernel <- kernel
res$deg <- deg
res$xeval <- xeval
## get info from formula
res$mf <- model.frame(formula,data)
datCla <- attr(attr(res$mf, "terms"),"dataClasses")
varNames <- names(datCla)[datCla=="numeric"]
stopifnot(length(varNames)==2)
res$Y <- varNames[1]
res$X <- varNames[2]
## sort x's
xo <- order(res$mf[,res$X])
res$mf <- res$mf[xo,]
res$weig <- weig[xo]
## xeval
if( is.null(xeval) )
res$xeval <- seq(min(res$mf[,res$X]),max(res$mf[,res$X]),len=xevalLen)
else
res$xeval <- sort(xeval)
if( is.null(res$bw) )
res$bw <- regCVBwSelC(data[,res$X],data[,res$Y],res$deg,
res$kernel,res$weig)
## regression estimation
res$lpFit <- locPolSmootherC(res$mf[,res$X], res$mf[,res$Y], res$xeval,
res$bw, res$deg, res$kernel, DET = TRUE, res$weig )
names(res$lpFit)[] <- c(res$X,res$Y,paste(res$Y,1:deg,sep=""),"xDen")
res$lpFit$xDen <- res$lpFit$xDen^(1/(deg+1))/(nrow(data)*res$bw)
## CI comp.(##Should depned on nu, to be able to)
nu <- 0
res$CIwidth <- computeRK(equivKernel(kernel,nu,deg),
lower=dom(res$kernel)[[1]], upper=dom(res$kernel)[[2]],
subdivisions = 25) * factorial(nu)^2
res$CIwidth <- res$CIwidth / ( nrow(data)*res$bw )
## residuals
res$residuals <- res$mf[,res$Y]-locLinSmootherC(res$mf[,res$X],
res$mf[,res$Y], res$mf[,res$X],res$bw, res$kernel,
res$weig)$beta0
## variance estimation
res$lpFit$var <- locCteSmootherC(res$mf[,res$X], res$residuals^2,
res$xeval,1.2*res$bw, res$kernel, res$weig)$beta0
## setupclass
class(res) <- "locpol"
return(res)
}
residuals.locpol <- function(object,...)
##
{
return( object$residuals )
}
fitted.locpol <- function(object,deg=0,...)
##
{
stopifnot(object$deg>=deg)
return( object$lpFit[,2+deg] )
}
summary.locpol <- function(object,...)
##
{
cat("\nKernel = \n\t")
print( body(object$kernel) )
cat("\n")
print( data.frame(n=nrow(object$mf), deg=object$deg, bw=object$bw,
ase=mean(resid(object)^2), row.names=" ") )
}
print.locpol <- function(x,...)
##
{
summary.locpol(x)
}
confInterval <- function(x)
##
{
plot(x$mf[,x$X],x$mf[,x$Y],pch="+",main="95% Conf. Int. for x-Points",
xlab=x$X,ylab=x$Y )
points(x$lpFit[,x$X],x$lpFit[,x$Y],type="l",col="darkgreen")
dev <- sqrt(x$CIwidth * x$lpFit$var/x$lpFit$xDen)
points(x$lpFit[,x$X],x$lpFit[,x$Y]+2*dev,type="l",col="green")
points(x$lpFit[,x$X],x$lpFit[,x$Y]-2*dev,type="l",col="green")
}
plot.locpol <- function(x,...)
##
{
par(ask=TRUE,cex=.6)
plot(x$mf[,x$X],x$mf[,x$Y],pch="+",main="Data and Regres.",
xlab=x$X,ylab=x$Y )
points(x$lpFit[,x$X],x$lpFit[,x$Y],type="l",col="blue")
ym <- max(x$lpFit$xDen)
plot(x$lpFit[,c(x$X,"xDen")], main="X dens.",type="l",ylim=c(0,1.25*ym),
xlab=x$X,ylab="den" )
ym <- max(x$lpFit$var)
plot(x$lpFit[,c(x$X,"var")], main="Var.",type="l",ylim=c(0,1.25*ym),
xlab=x$X,ylab="var" )
confInterval(x)
par(ask=FALSE)
}
# mvNoparEst.R
# Multivariate density and regression.
# NOTAS:
# . Check normalization constant !!
# ERRORS:
##
## some bivariate kernels,
##i
epaK2d <- function(x) (2/pi) * (1-x[[1]]^2 -x[[2]]^2)*( (x[[1]]^2+x[[2]]^2)< 1 )
gauK2d <- function(x) (1/(2*pi)) * exp(-0.5*(x[[1]]^2 + x[[2]]^2) ) ##check this, not sure !!!
##
## ... a bandwidth selector ??
##
##myDR <- function(x) diff(range(x))
mayBeBwSel <- function(X,prop=.45)
##
# X = data(data.frame or matrix), only first two cols are used.
{
## uses 'prop' of the range of points.
return( diag(prop*c(diff(range(X[,1])), diff(range(X[,2]))) ))
}
##
## density estimation
##
bivDens <- function(X,weig=rep(1,nrow(X)),K=epaK2d,H=mayBeBwSel(X) )
##
# X = data(data.frame or matrix), only first two cols are used.
# weig = weigths for each point.
# K = bivariate kernel.
# H = bandwidht matrix, by default
{
stopifnot( length(weig)==nrow(X) )
weig <- weig/sum(weig)
## compute kernel dens. est. for each data point.
denEst <- function(xa)
{
we <- apply(X[,1:2],1,function(x)K( (x-xa) %*% solve(H)))
return( sum(we*weig) )
}
## apply it to each data point
res <- list( X=X[,1:2], H=H, estFun=denEst )
class(res) <- "bivNpEst"
return(res)
}
##
## Regression
##
bivReg <- function(X,Y,weig=rep(1,nrow(X)),K=epaK2d,H=mayBeBwSel(X) )
##
# X = X data(data.frame or matrix), only first two cols are used.
# Y = Y data,a vector.
# weig = weigths for each point.
# K = bivariate kernel.
# H = bandwidht matrix.
{
stopifnot( length(weig)==nrow(X), length(Y)==nrow(X))
weig <- weig/sum(weig)
## compute NW est. for each data point.
regEst <- function(xa)
{
we <- apply(X[,1:2],1,function(x)K( (x-xa) %*% solve(H)))
return( sum(we*weig*Y)/sum(we*weig) )
}
## apply it to each data point
res <- list( X=X[,1:2], Y=Y, H=H, estFun=regEst )
class(res) <- "bivNpEst"
return(res)
}
predict.bivNpEst <- function(object,newdata=NULL,...)
##
# object = "bivNpEst" object
# newdata = new data where the density should be predicted
{
if( is.null(newdata) )
newdata <- object$X
else{
stopifnot( is.data.frame(newdata),all(names(object$X) %in% names(newdata)) )
}
return( apply(newdata[,names(object$X)],1,function(xa) object$estFun(xa)) )
}
plotBivNpEstOpts <- list(
pathLen=10, phi=30, theta=15, r=sqrt(3),
ticktype="detailed", nticks=3, main="NP est."
)
plot.bivNpEst <- function(x,...)
##
# x = a 'bivNpEst' object
# ... = 'persp' options plus 'plRng' and 'pathLen':
# plRng = matrix with the range for x and y in each row.
# pathLen = number of points in the grid
{
## default options
mc <- match.call()
parIdx <- na.omit( match(names(mc)[-1],names(plotBivNpEstOpts)) )
if( length(parIdx)>0 )
plotBivNpEstOpts <- plotBivNpEstOpts[-parIdx]
for( i in names(plotBivNpEstOpts) )
mc[[i]] <- plotBivNpEstOpts[[i]]
if( is.null(mc$plRng) )
plRng <- rbind(range(x$X[,1]),range(x$X[,2]))
else{
plRng <- eval(mc$plRng)
mc$plRng <- NULL
}
pathLen <- eval(mc$pathLen)
mc$pathLen <- NULL
if( is.null(mc$xlab) )
mc$xlab <- names(x$X)[1]
if( is.null(mc$ylab) )
mc$ylab <- names(x$X)[2]
if( is.null(mc$zlab) )
mc$zlab <- ifelse( is.null(x$Y),"den","reg")
##plot
mc[[1]] <- persp
mc[["x"]] <- seq(from=plRng[[1,1]],to=plRng[[1,2]],len=pathLen)
mc[["y"]]<- seq(from=plRng[[2,1]],to=plRng[[2,2]],len=pathLen)
dd <- expand.grid(mc[["x"]],mc[["y"]])
names(dd) <- names(x$X)
mc[["z"]]<- predict(x, dd)
dd$den <- mc[["z"]]
mc[["z"]] <- matrix(mc[["z"]],nrow=pathLen)
eval(mc)
return(dd)
}
#
# thumbBw.R
# Rule of thumb bandwith selector for loc. pol.
# NOTES:
# Pesos no implementados...
# ERRORS:
##source('kernels.R')
compDerEst <- function(x,y,p,weig=rep(1,length(y)))
##OK
# x = x data.
# y = y data.
# p = loc. pol. degree.
# weig = weigths used in the estimation
{
## compute global pol. regression
xnam <- paste("x^", 2:(p+3), sep="")
xnam <- paste("I(",xnam, ")")
fmla <- as.formula(paste("y ~ 1 + x + ", paste(xnam, collapse= "+")))
lmFit <- lm(fmla,data=data.frame(x,y),weights=weig)
## compute (p+1)-der
cp <- cumprod(1:(p+3))
coef <- coefficients(lmFit)
der <- coef[[p+2]]*cp[[p+1]]+coef[[p+3]]*cp[[p+2]]*x+
coef[[p+4]]*cp[[p+3]]*x^2/2
##der <-
res <- data.frame(x,y,res=residuals(lmFit)*weig,der)
return( res )
}
thumbBw <- function(x,y,deg,kernel,weig=rep(1,length(y)))
##OK
# x = x data.
# y = y data.
# deg = loc. pol. degree.
# kernel = kenrel used in estimation.
# weig = weigths used in the estimation
{
k <- 3
## compute reg. residuals and derivatives
rd <- compDerEst(x,y,deg,weig)
denom <- sum( rd$der^2 )
numer <- mean( rd$res^2 )
##constante kernel, eta, p
cte <- cteNuK(0,deg,kernel,lower=dom(kernel)[[1]],
upper=dom(kernel)[[2]],subdivisions=100)
##valor de h
res <- cte * (numer/denom)^(1/(2*deg+k))
return( res )
}
pluginBw <- function(x,y,deg,kernel,weig=rep(1,length(y)))
##OK
# x = x data.
# y = y data.
# deg = loc. pol. degree.
# kernel = kenrel used in estimation.
# weig = weigths used in the estimation
# Only valid for odd degrees, base on Fand&Gijbels book.
{
stopifnot(deg%%2==1)
## thumb bandwithd
thBw <- thumbBw(x, y, deg, kernel)
## compute loo res.
regComp <- looLocPolSmootherC(x, y, thBw, deg, kernel, weig)
res <- (y-regComp$beta0)*weig
numer <- mean( res^2 )
## compute loo der.
thBwBis <- thumbBw(x, y, deg+2, kernel)
derBw <- adjNuK(deg+1,deg+2,kernel,lower=dom(kernel)[[1]],
upper=dom(kernel)[[2]],subdivisions=50) * thBwBis
regCompBis <- looLocPolSmootherC(x, y, derBw, deg+2, kernel, weig)
cp <- cumprod(1:(deg+2))
der <- regCompBis[,1+deg+2]/cp[[deg+1]]
denom <- sum( der^2 )
##constante kernel, eta, p
cte <- cteNuK(0,deg,kernel,subdivisions=100)
##valor de h
res <- cte * (numer/denom)^(1/(2*deg+3))
return(res)
}
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.