L = 100
X1 = seq(9,21,length=L)
X2 = seq(3,17,length=L)
ev <- function(x1,x2){ return(90 + 150*x1 + 105*x2 - 4*x1^2 - 3*x2^2 - 3*x1*x2) } Y = matrix(NA,ncol=L,nrow=L) for(i in 1:L){ for(j in 1:L){ Y[i,j] = ev(X1[i],X2[j]) } } contour(X1,X2,Y) persp(X1,X2,Y) # Generate noisy observations: y = Y + matrix(10*rnorm(L*L),ncol=L,nrow=L) + matrix(50*runif(L*L,-1,1),ncol=L,nrow=L) contour(X1,X2,Y) contour(X1,X2,y,add=TRUE,col=2) persp(X1,X2,y) # polynomial regression: gg = expand.grid(X2,X1) dat = data.frame(y=c(y), x1=gg[,2], x2=gg[,1]) loo = loess(y~., data=dat) yhat = matrix(loo$fitted,ncol=L,nrow=L) contour(X1,X2,Y) contour(X1,X2,y,add=TRUE,col=2) contour(X1,X2,yhat,add=TRUE,col=4) persp(X1,X2,yhat) persp(X1,X2,y-yhat) # model response as a function of theta: model <- function(th,x1,x2){ return(th[1] + th[2]*x1 + th[3]*x2 - th[4]*x1^2 - th[5]*x2^2 - th[6]*x1*x2) }