# Tests using the functions in gbasis par(mfrow=c(2,4)) # Generate an example data set. set.seed(2) x.grid <- seq(0,1,0.005) t.grid <- sin(1+x.grid^2) xa <- runif(1050) ta <- sin(1+xa^2) + rnorm(1050,0,0.03) xt <- xa[1:50] tt <- ta[1:50] # Try out regularized estimation with penalty set by cross-validation. try.lambda <- c(0.001,0.01,0.1,1,10) try.s <- c(0.02,0.1,0.5,2.5) cat("\nCross-validation array:\n\n") print (round(sqrt(cross.val.array(xt,tt,try.s,try.lambda,5)),4)) cat("\nRoot mean square error on test cases:\n\n") print (round(sqrt(val.array(xa,ta,try.s,try.lambda,51:1050)),4)) for (s in c(0.1,0.5)) { for (lambda in c(0.01,1)) { plot (xt, tt, pch=20, ylab="", xlab=paste("Gaussian basis function fit, s =",s," lambda =",lambda)) lines(x.grid, t.grid, lwd=2, col="gray") w <- penalized.least.squares (Phi(xt,s), tt, lambda) p.grid <- predictions (Phi(x.grid,s), w) lines(x.grid ,p.grid, lwd=2, col="black") } } # Try out Bayesian inference with hyperparameters from marginal likelihood. try.sigma <- seq(0.020,0.045,by=0.001) omega0 <- 10 for (s in c(0.1,0.5)) { cat("\nLog marginal likelihoods, s = ",s,":\n\n") try.omega <- if (s==0.1) seq(0.020,0.070,by=0.002) else seq(0.050,0.550,by=0.020) Phi.s <- Phi(xt,s) ml <- matrix(NA,length(try.sigma),length(try.omega)) rownames(ml) <- paste("sigma=",try.sigma,sep="") colnames(ml) <- paste("omega=",try.omega,sep="") best.ml <- -Inf for (si in 1:length(try.sigma)) { for (oi in 1:length(try.omega)) { ml[si,oi] <- marg.like (Phi.s, tt, try.sigma[si], try.omega[oi], omega0) if (ml[si,oi]>best.ml) { best.sigma <- try.sigma[si] best.omega <- try.omega[oi] best.ml <- ml[si,oi] } } } print(round(ml)) cat("\nBest log marginal likelihood:",best.ml, " best sigma:",best.sigma," best omega",best.omega, " corresponding lambda:",best.sigma^2/best.omega^2,"\n") contour(try.sigma,try.omega,exp(ml-best.ml),xlab="sigma",ylab="omega") title(paste("s =",s)) plot (xt, tt, pch=20, ylab="", xlab=paste("Bayesian fit, s =",s," sigma =",best.sigma, " omega =",best.omega)) lines(x.grid, t.grid, lwd=2, col="gray") w <- posterior.mean (Phi(xt,s), tt, best.sigma, best.omega, omega0) p.grid <- predictions (Phi(x.grid,s), w) lines(x.grid ,p.grid, lwd=2, col="black") }