par(mai=c(0,0,0,0)) # error standard deviation: sigma <- 0.1 # true regression function, put in dataframe fframe: f <- function(x, y) 5*x*(x-1)*y*(y-1)*sin(11*x)*sin(11*y) #f <- function(x, y) 5*x*(x-1)*y*(y-1)*sin(8*x)*sin(8*y) v <- seq(0,1, length=51) w <- seq(0,1, length=51) fframe <- expand.grid(v=v, w=w) fframe\$f <- f(fframe\$v, fframe\$w) persp(v,w, matrix(fframe\$f, 51, 51), theta = 30, phi = 30, expand = 0.5, zlab="", xlab="x", ylab="y", main="") -> res maxindex <- which(fframe\$f==max(fframe\$f)) xfmax <- fframe\$v[maxindex] yfmax <- fframe\$w[maxindex] #points(trans3d(xfmax,yfmax, f(xfmax, yfmax), pmat=res), col=2) # initial, uniformly spaced data, put in dataframe obsframe: # totaal n1*n1 observaties n1 <- 25 x <- seq(0,1, length=n1) y <- seq(0,1, length=n1) obsframe <- expand.grid(x=x, y=y) obsframe\$z <- f(obsframe\$x, obsframe\$y) + sigma*rnorm(n1*n1) #persp(x,y, matrix(obsframe\$z, n1, n1), theta = 30, phi = 30, expand = 0.5, zlab="z", xlab="x", ylab="y", main="noisy observations") maxindex <- which(fframe\$f==max(fframe\$f)) xfmax <- fframe\$v[maxindex] yfmax <- fframe\$w[maxindex] points(trans3d(obsframe\$x,obsframe\$y, obsframe\$z, pmat=res), col="gray") # local linear regression to find stage one estimator: m <- loess(obsframe\$z ~obsframe\$x + obsframe\$y, degree =1, span=.0855) fit <- predict(m) persp(x,y, matrix(fit, n1, n1), theta = 30, phi = 30, expand = 0.5, zlab="", xlab="x", ylab="y", main="") maxindex <- which(fit==max(fit)) xmutilde <- obsframe\$x[maxindex] ymutilde <- obsframe\$y[maxindex] #points(trans3d(xmutilde,ymutilde, f(xfmax, yfmax), pmat=res), col=2) # new data n3 <- 70 delta <- .1 # new designpoints: x0 <- xmutilde y0 <- ymutilde x1 <- x0+delta y1 <- y0 x2 <- x0 y2 <- y0+delta x3 <- x0-delta y3 <- y0 x4 <- x0 y4 <- y0-delta x5 <- x0+delta y5 <- y0-delta x6 <- x0+delta y6 <- y0+delta x7 <- x0-delta y7 <- y0+delta x8 <- x0-delta y8 <- y0-delta # regressors: p <- rep(c(x0, x1, x2, x3, x4, x5, x6, x7, x8), n3) q <- rep(c(y0, y1, y2, y3, y4, y5, y6, y7, y8), n3) r <- rep(c((x0)^2, (x1)^2, (x2)^2, (x3)^2, (x4)^2, (x5)^2, (x6)^2, (x7)^2, (x8)^2), n3) s <- rep(c((y0)^2, (y1)^2, (y2)^2, (y3)^2, (y4)^2, (y5)^2, (y6)^2, (y7)^2, (y8)^2), n3) t <- rep(c(x0*y0, x1*y1, x2*y2, x3*y3, x4*y4, x5*y5, x6*y6, x7*y7, x8*y8), n3) # new observations: newz <- f(p,q) + sigma*rnorm(9*n3) points(trans3d(c(x0, x1, x2, x3, x4, x5, x6, x7, x8), c(y0, y1, y2, y3, y4, y5, y6, y7, y8), newz, pmat=res), col=8) points(trans3d(xmutilde,ymutilde, f(xfmax, yfmax), pmat=res), col=2, type="b",pch=19) a <- lm(newz ~ p+q+r+s+t)\$coefficients # find muhat, the argmax of the quadratic surface A <- matrix(data=c(2*a[4], a[6], a[6], 2*a[5]), nrow=2, ncol=2, byrow=TRUE) b <- -c(a[2], a[3]) u <- solve(A,b) xmuhat <- as.numeric(u[1]) ymuhat <- as.numeric(u[2]) # plot of quadratic surface: g <- function(k,l) a[1] + a[2]*k + a[3]*l + a[4]*k^2+a[5]*l^2+a[6]*k*l gframe <- expand.grid(v=v, w=w) gframe\$g <- g(gframe\$v, gframe\$w) persp(v,w, matrix(gframe\$g, 51, 51), theta = 30, phi = 30, expand = 0.5, zlab="", xlab="x", ylab="y", main="") -> res # did we make things worse? problem <- (abs(xmutilde - xmuhat) > delta) | (abs(ymutilde - ymuhat) > delta) if (problem) { xmuhat <- xmutilde ymuhat <- ymutilde } # plot the final estimator: points(trans3d(xmuhat,ymuhat, g(xmuhat, ymuhat), pmat=res), col="green", pch=19)