#Price of stamps in cents
y = c(.05,.06,.08,.10,.13,.15,.20,.22,.25,.29,.32,.33,.34,.37,.39,.41,.42,.44,.45,.46)*100

#years since 1960
x = c(63,68,71,74,75,78,81,85,88,91,95,99,101,102,106,107,108,109,112,113)-60

plot(x,y)
fit = lm(y~x)
print(fit)
print(fit$coefficients)

# by hand
xb = mean(x)
yb = mean(y)
sxx = mean((x-xb)^2)
sxy = mean((x-xb)*(y-yb))
b1 = sxy/sxx
b0 = yb - b1*xb


# sample fitting code ----
x <- 1:10
y <- .25*(x*(x-10)+27) + c(-.7,.7)
par('mar'=c(4,4,.5,.5))
plot(x,y, xlim=c(0,11), ylim=c(-2,13),lwd=2,cex.lab=1.5)
fit1 <- lm( y~x )
fit2 <- lm( y~poly(x,2) )
fit3 <- lm( y~poly(x,9) )
xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1, data.frame(x=xx)), col='blue',lwd = 2)
lines(xx, predict(fit2, data.frame(x=xx)), col='maroon',lwd=2)
lines(xx, predict(fit3, data.frame(x=xx)), col='orange',lwd=2)
#compute R^2
yhat = predict(fit1,data.frame(x=x))
R1sq = cor(y,yhat)
yhat = predict(fit2,data.frame(x=x))
R2sq = cor(y,yhat)
yhat = predict(fit3,data.frame(x=x))
R3sq = cor(y,yhat)
print(c(R1sq,R2sq,R3sq))
# MORE FITS
#fit4 <- lm( y~offset(x) -1 )
#library(splines)
#fit5 <- lm( y~ns(x, 3) )
#fit6 <- lm( y~ns(x, 9) )
#fit7 <- lm( y ~ x + cos(x*pi) )
#lines(xx, predict(fit4, data.frame(x=xx)), col='purple')
#lines(xx, predict(fit5, data.frame(x=xx)), col='orange')
#lines(xx, predict(fit6, data.frame(x=xx)), col='grey')
#lines(xx, predict(fit7, data.frame(x=xx)), col='black')