##
## ST4060 / ST6040
## Review script – Lecture 15
## (complementing grid search…)
## Eric Wolsztynski
##
while( dev.next()>1 ){ dev.off() }
rm(list=ls(pos=.GlobalEnv), pos=.GlobalEnv)
# airquality
dat = na.omit(airquality)
x = dat$Temp
y = dat$Ozone
par(font=2, font.lab=2, font.axis=2) # use bold fonts in figure
plot(x, y, pch=20, xlab=”Temperature (degrees F)”, ylab=”Ozone (ppb)”,
xlim=c(55,105))
nl.crit <- function(x, y, a=1, b=0.01){ return( mean( (y - a*exp(b*x))^2 ) ) } # GRID SEARCH: Pa = Pb = 100 a.grid = seq(0,2,length=Pa) b.grid = seq(0,1,length=Pb) # try this one first... # try again with: # a.grid = seq(-0.75,1.2,length=Pa) # b.grid = seq(0,.05,length=Pb) # then this one... # b.grid = seq(.04,.07,length=Pb) # try last crit.values = matrix(NA, nrow=Pa, ncol=Pb) for(i in 1:Pa){ for(j in 1:Pb){ crit.values[i,j] = nl.crit(x,y,a.grid[i],b.grid[j]) } } ij.best = which.min(crit.values) ij.best = arrayInd(ij.best, .dim=dim(crit.values)) # grid search solutions: a.best = a.grid[ij.best[1]] b.best = b.grid[ij.best[2]] fitted = a.best * exp(b.best*x) par(font=2, font.lab=2, font.axis=2) # use bold fonts in figure par(mfrow=c(2,2)) # 1st plot plot(x, y, pch=20, xlab="Temperature (degrees F)", ylab="Ozone (ppb)", xlim=c(55,105)) points(x, fitted, col='orange', pch=20) # 2nd plot plot(a.grid, crit.values[,ij.best[2]], t='b') # 3rd plot plot(b.grid, crit.values[ij.best[1],], t='b') # 4th plot persp(a.grid, b.grid, crit.values, theta=30, phi=10) dev.new() persp(a.grid, b.grid, crit.values, theta=-180, phi=0) contour(a.grid, b.grid, crit.values, nlevels=100) nlso = nls(y~a*exp(b*x), start=list(a=a.best, b=b.best)) a.best b.best plot(x, y, pch=20, xlab="Temperature (degrees F)", ylab="Ozone (ppb)", xlim=c(55,105)) points(x, fitted(nlso), pch=20, col='orange') # try again with: # a.grid = seq(-0.75,1.2,length=Pa) # b.grid = seq(0,.05,length=Pb) # then this one... # b.grid = seq(.04,.07,length=Pb) # try last