rm(list = ls()) # Simulated responses from # y = 5 + cos((2*w1 + w2 - 15)/10) + cos((w1 - 2*w2 + 5)/10) + .15*z # where z is standard normal. ################################################################## W1 = seq(0,20,by=.5) #sequences W1,W2 with elements w1,w2 W2 = W1 meany = function(w1, w2) 5 + cos((2*w1 + w2 - 15)/10) + cos((w1 - 2*w2 + 5)/10) maty = outer(W1, W2, function(w1, w2) meany(w1,w2)) persp(W1,W2,maty,ticktype="detailed",col = "red") windows() contour(W1,W2,maty) centre = c(10,10) iter = 0 text(centre[1], centre[2], iter) p.inter = 1 p.curv = 1 OUT = rbind(c(centre, p.inter, p.curv)) colnames(OUT) = c("w1", "w2", "p.inter", "p.curv") ############ Run until curvature is significant ################# #Code w1 and w2: W1 = function(x1) 2*x1 + centre[1] W2 = function(x2) 2*x2 + centre[2] # First order design: 2^2 factorial and 5 centre points x1 = c(-1,-1,1,1,0,0,0,0,0) x2 = c(-1,1,-1,1,0,0,0,0,0) w1 = W1(x1) w2 = W2(x2) y = vector(length = length(x1)) for(i in 1:length(y)) y[i] = meany(w1[i],w2[i]) + .15*rnorm(1) mat1 = cbind(x1,x2,w1,w2,y) mat1 y = mat1[,5] x = cbind(x1, x2, x1*x2, x1*x1) colnames(x) = c("x1", "x2", "inter'n", "curv") fit1 = lsfit(x,y) out = ls.print(fit1) p.inter = out$"coef.table"[[1]][4,4] p.curv = out$"coef.table"[[1]][5,4] nf = 4 nc = 5 yf = y[1:nf] yc = y[(nf+1):(nf+nc)] sigmasqd = var(yc) t.curv = (mean(yf) - mean(yc))/sqrt(sigmasqd*(1/nf + 1/nc)) p.curv = 2*pt(abs(t.curv), nc-1, lower = FALSE) p.curv beta = fit1$coef[2:3] beta = beta/sqrt(sum(beta^2)) beta # Make observations at (x1,x2) = lambda*beta lambda = 1:10 #lambda = lambda/10 ## Use this if a loop yields no change becuase lambda = 1 is too large x1 = lambda*beta[1] x2 = lambda*beta[2] w1 = W1(x1) w2 = W2(x2) y = vector(length = length(x1)) for(i in 1:length(lambda)) y[i] = meany(w1[i],w2[i]) + .15*rnorm(1) mat2 = cbind(lambda,x1,x2,w1,w2,y) mat2 = rbind(c(0,0,0,centre,mean(yc)), mat2) colnames(mat2) = c("lambda", "x1", "x2", "w1", "w2", "y") mat2 # Compute new centre: ally = mat2[,6] bestrow = which(ally == max(ally)) centre = mat2[bestrow,4:5] iter = iter + 1 text(centre[1], centre[2], iter) OUT = rbind(OUT, c(centre, p.inter, p.curv)) out round(OUT,3) ################################################################## ################################################################## meany = function(w1, w2) 5 + cos((2*w1 + w2 - 15)/10) + cos((w1 - 2*w2 + 5)/10) #centre = c(5.963339, 5.567325 ) W1 = function(x1) 2*x1 + centre[1] W2 = function(x2) 2*x2 + centre[2] # Second order design: 2^2 factorial and 5 centre points and 4 axial points alpha = sqrt(2) x1 = c(-1,-1,1,1,0,0,0,0,0, alpha, -alpha,0,0) x2 = c(-1,1,-1,1,0,0,0,0,0,0,0,alpha, -alpha) # Note that 9 of these runs have already been made and we could have # (but don't here) re-use these results. w1 = W1(x1) w2 = W2(x2) y = vector(length = length(x1)) for(i in 1:length(y)) y[i] = meany(w1[i],w2[i]) + .15*rnorm(1) mat3 = cbind(x1,x2,w1,w2,y) mat3 x = cbind(x1, x2, x1*x2, x1*x1, x2*x2) colnames(x) = c("x1", "x2", "inter'n", "quad x1", "quad x2") fit3 = lsfit(x,y) out = ls.print(fit3) out coef = fit3$coef[-1] b = c(coef[1], coef[2]) B11 = coef[4] B22 = coef[5] B12 = .5*coef[3] B = rbind(c(B11, B12), c(B12, B22)) x.stat = -.5*solve(B,b) w1 = W1(x.stat[1]) w2 = W2(x.stat[2]) cat("Estimated stationary point is w1 = ", round(w1,3), "and w2 = ", round(w2,3), "\n") cat ("Eigenvalues of Hessian are", eigen(2*B, only.values = T)$values, "\n")