x <- seq(-1,1,,50)
y <- (f.true <- pnorm(2*x)) + rnorm(50)/10
## specify pointwise constraints (boundary conditions)
con <- rbind(c( 1,min(x),0), # f(min(x)) >= 0
c(-1,max(x),1), # f(max(x)) <= 1
c(0, 0, 0.5))# f(0) = 0.5
## obtain the median regression B-spline using automatically selected knots
Rbs <- cobs(x,y,constraint="increase",pointwise=con)
Rbs
plot(x,y)
lines(predict(Rbs), col = 2, lwd = 1.5)
lines(spline(x,f.true), col = "gray40")
Rbsub <- cobs(x,y,constraint="increase",pointwise=con, n.sub = 45)summary(Rbsub)
lines(predict(Rbsub), col = 4, lwd = 1)
## compute the median smoothing B-spline using automatically chosen lambda
Sbs <- cobs(x,y,constraint="increase",pointwise=con,lambda=-1)
Sbs
plot(Sbs$pp.lambda[-1], Sbs$sic[-1], log = "x",
main = "SIC ~ lambda", xlab = expression(lambda), ylab = "SIC")
axis(1, at = Sbs$lambda, label = expression(hat(lambda)),
col.axis = 2, mgp = c(3, 0.5, 0))
Sb1 <- cobs(x,y,constraint="increase",pointwise=con,lambda=-1, degree=1)
summary(Sb1)
pxx <- predict(Sb1, xx <- seq(-1.2, 1.2, len = 201), interval = "both")
plot(x,y, main = deparse(Sb1$call),
xlim = range(xx), ylim = range(y, pxx[,"fit"]))
lines(pxx, col = 2)
rug(Sb1$knots, col = 4, lwd = 1.6)# (too many knots)
matlines(pxx[,1], pxx[,-(1:2)], col= rep(c("green","blue"),c(2,2)), lty=2)
Run the code above in your browser using DataLab