################################################################################# ## Spline smoothing illustrated as BLP prediction ################################################################################# source("/net/roxbee/home/pmcc/public_html/courses/regress.R") nobs <- 100 x1 <- 10*runif(nobs) # observation points x0 <- seq(-1, 15, 0.1) # prediction points (ordered for plotting) npred <- length(x0) x <- c(x0, x1) obs <- c(rep(FALSE, npred), rep(TRUE, nobs)) mu <- log(2+x) + 0.5* sin(3*log(2+x)) # arbitrary non-linear `random' function sigma <- 0.3 nn <- length(x) y <- mu + rnorm(nn) * sigma d <- abs(outer(x, x, "-")) # distance matrix lambda <- 10 K <- exp(-d/lambda) # exponential covariance function fit0 <- regress(y[obs]~x[obs]) fit1 <- regress(y[obs]~x[obs], ~K[obs, obs]) 2*(fit1$llik - fit0$llik) summary(fit1) fitted <- fit1$beta[1] + fit1$beta[2]*x # (fitted[obs] = fit1$fitted) blp <- fitted[!obs] + fit1$sigma[2]*K[!obs, obs] %*% fit1$W %*% (y[obs] - fitted[obs]) plot(x, y, cex=(1+obs)/3) lines(x0, blp, col="blue") # best linear predictor lines(x0, mu[!obs], col="green", lty=2) # target mean sum((y[!obs] - blp)^2)/npred # mean squared prediction error ################################################################################## #### repeat with different K ################################################################################## K <- d^2*log(d+diag(nobs+npred))/100 K <- d^3/10^3 fit0 <- regress(y[obs]~x[obs]) fit1 <- regress(y[obs]~x[obs], ~K[obs, obs]) 2*(fit1$llik - fit0$llik) fitted <- fit1$beta[1] + fit1$beta[2]*x blp <- fitted[!obs] + fit1$sigma[2]*K[!obs, obs] %*% fit1$W %*% (y[obs] - fitted[obs]) plot(x, y, cex=(1+obs)/3) lines(x0, blp, col="blue") # best linear predictor lines(x0, mu[!obs], col="green", lty=2) # target mean sum((y[!obs] - blp)^2)/npred # mean squared prediction error ################################################################################## K1 <- -d # for linear splines K2 <- d^2*log(d+diag(n))/10^2 # for quadratic splines K3 <- d^3/10^3 # for cubic splines