## ——————————————
##
## ST4060 / ST6015 / ST6040
## R script – Wedn 27 Oct 2021 lecture
## Eric Wolsztynski
##
## ——————————————
rm(list=ls())
attach(pressure) # NB: you should usually not use attach() 😉
x = temperature
y = pressure
plot(x, y, pch=20, cex=2)
# We can use a linear estimation procedure (lm() here)
# to fit the polynomial, because it is linear in the parameters
# of interest:
lmo = lm(y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5))
lines(x, fitted(lmo), lwd=4, col=4)
# The polynomial does an increasingly better job at fitting
# as we increase its degree… which leads to “overfitting”
# (sometimes that isn’t a problem though!).
# How to fit a nonlinear model using nls():
nlreg <- nls(y~exp(a+b*x), start=list(a=0,b=.02))
summary(nlreg)
cc = coef(nlreg)
# plot the pattern:
par(font=2, font.lab=2, font.axis=2)
plot(x, y, pch=20, cex=2,
xlab="Temperature (Celsius)",
ylab="Vapor pressure (ml of mercury)",
main="Example: polynomial regression")
curve(exp(cc[1]+cc[2]*x), col='blue', add=T, lwd=3)
# How to fit any model using optim():
crit <- function(theta,x,y){
# Must have theta as 1st parameter and return a single value...
# Here we impolement a simple linear regression model
# for the sake fo the example:
return( sum( (y-x*theta)^2 ) )
}
thbar = 1.8
x = rep(c(1,3,7,8), len=100)
y = x*thbar + rnorm(x)
(optim.out = optim(par=c(1), fn=crit, x=x, y=y, method='L-BFGS-B',
lower=c(0.01), upper=c(3.05)))
# Try optim() with a nonlinear moedl (e.g. the exponential model)!