# demo/constrained.R In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

## $Id: lc_k1_mean.R,v 1.11 2013/12/01 16:07:16 jracine Exp jracine$

## We illustrate constrained kernel estimation using the local
## constant estimator where the constraints are l(x) <= \hat g(x) <=
## u(x) (see Du, P. and C. Parmeter and J.S. Racine (2013),
## "Nonparametric Kernel Regression with Multiple Predictors and
## Multiple Shape Constraints," Statistica Sinica, Volume 23, Number
## 3, 1343-1372).

rm(list=ls())

## Load required packages, set options.

library(np)
library(crs)
options(np.tree=TRUE,np.messages=FALSE)

## Set the kernel function.

ckertype <- "epanechnikov"

## Set the order of the local polynomial p' (p=0 is the
## Nadaraya-Watson estimator, p=1 local linear, p=2 local quadratic
## etc.). If p < 0 both p and the bandwidth (bw) are determined by
## cross-validation, otherwise only bw is determined.

p <- 1

## Simulate a sample of data. We can control signal/noise ratio by
## multiplying epsilon by sd(dgp), so rnorm(n,sd=...) can be set to
## sd=(.25,.5,1,2) which would yield an R-squared for the Oracle model
## of (.95,.8,.5,and .2). Changing n changes the sample size,
## modifying cos() the data generating process (dgp).

n <- 1000
x <- sort(runif(n))
dgp <- cos(2*pi*x)
y <- dgp + sd(dgp)*rnorm(n,sd=0.25)

## Set up bounds for the quadratic program. We are going to require
## the lower and upper constraints l(x) and u(x) for g(x). Here the
## upper bound is a constant, the lower bound depends on x.

lower <- -1.5+2*x
upper <- rep(0.5,n)
if(any(lower>upper)) stop(" constraints are inconsistent")

## X (data frame of regressors) and y are passed below, so if you add
## extra regressors simply add them to X here and to the formula (used
## to compute cross-validated smoothing parameters) and be done.

X <- data.frame(x)

formula.glp <- formula(y~x)

model.glp <- crs:::npglpreg(formula=formula.glp,
cv=ifelse(p>=0,"bandwidth","degree-bandwidth"),
degree=ifelse(p>=0,rep(p,NCOL(X)),rep(0,NCOL(X))),
ckertype=ckertype,
nmulti=min(NCOL(X),5))

bw <- model.glp$bws p <- model.glp$degree

## W is the local polynomial model matrix (includes an intercept).

W <- crs:::W.glp(xdat=X,
degree=rep(p,NCOL(X)))

## Generate the matrix of kernel weights using data-driven bandwidths
## that are optimal for the unconstrained model.

K <- npksum(txdat=X,
bws=bw,
ckertype=ckertype,
return.kernel.weights=TRUE)$kw ## Compute the matrix for which \hat g(x|p) = t(A)%*%p (i.e. the local ## polynomial fit). A <- sapply(1:n,function(i){W[i,,drop=FALSE]%*%chol2inv(chol(t(W)%*%(K[,i]*W)))%*%t(W)*y*K[,i]}) ## Compute the unconstrained estimator. p.u <- rep(1,n) fit.unres <- t(A)%*%p.u ## Solve the quadratic program. The function solve.QP in the quadprog ## package solves the problem min (p-p.u)'(p-p.u) subject to the ## constraints Amat^T p >= bvec. Note that we construct Amat to ## contain the constraints a) ) A^T p >= lower, and b) -A^Tp >= -upper ## (here A^T = \hat g(x|p)). output.QP <- solve.QP(Dmat=diag(n), dvec=rep(1,n), Amat=cbind(A,-A), bvec=c(lower,-upper)) if(is.nan(output.QP$value)) stop(" solve.QP was unable to find a solution: adjust constraints and restart")

p.hat <- output.QP\$solution

## Compute the constrained estimator.

fit.res <- t(A)%*%p.hat

## Plot the DGP, unrestricted, and restricted fits along with the
## constraints and data.

ylim <- c(min(dgp,lower,fit.unres,fit.res,y),
max(dgp,upper,fit.unres,fit.res,y))
subtext <- paste("Local polynomial order = ", p,", bandwidth = ",formatC(bw,digits=3,format="f"),", n = ",n,sep="")
plot(x,y,cex=.25,ylab="g(x)",ylim=ylim,col="grey",sub=subtext)
lines(x,dgp,col=1,lty=1)
lines(x,fit.unres,col=2,lty=1,lwd=2)
lines(x,fit.res,col=3,lty=1,lwd=2)
lines(x,lower,lty=1,col="lightgrey")
lines(x,upper,lty=1,col="lightgrey")

legend("topleft",
c("DGP","Unconstrained [g(x)]","Constrained [g(x|p)]","Constraints [l(x) <= g(x|p) <= u(x)]"),
lty=c(1,1,1,1),
col=c(1:3,"lightgrey"),
bty="n",
cex=0.75)
`

## Try the np package in your browser

Any scripts or data that you put into this service are public.

np documentation built on March 31, 2023, 9:41 p.m.