## Podiplomski studij statistike, 2002/2003
## V. Batagelj: Informacijska tehnologija v analizi podatkov
## Optimizacija v R-ju
## 10.jan 03 / 05.jan 03
##
# ------------------------------------------------------------------
## Linearna optimizacija

a <- read.table(file="playfair.dat",header=TRUE,sep="",row.names=1)
plot(a)
attach(a); plot(diameter,population)
b <- lm(population ~ diameter)
b
summary(b)
plot(b)
plot(diameter,population)
abline(b,col="red")
pb <- predict(b)
points(diameter,pb,pch=16,col="red")

plot(diameter,population)
c <- lm(population ~ 1 + diameter + I(diameter^2))
x <- seq(5,41,1)
pc <- function(x){coef(c)[3]*x^2 + coef(c)[2]*x + coef(c)[1]}
lines(spline(x,pc(x)),col="red")
points(diameter,pc(diameter),pch=16,col="red")

# ------------------------------------------------------------------
## Nelinearna optimizacija

f4 <- function(x){x^4 - 14*x^3 + 60*x^2 - 70*x}
curve(f4,0,7)
m4 <- optimize(f4,interval=c(5,7),tol=0.000001)
lines(c(m4$min,m4$min,-100),c(-100,m4$obj,m4$obj),col="red")

fr <- function(x) { ## Rosenbrock Banana function
  x1 <- x[1]; x2 <- x[2];
  100 * (x2 - x1^2)^2 + (1 - x1)^2
}
gr <- function(x) { ## Gradient of `fr'
  x1 <- x[1]; x2 <- x[2]
  c(-400*x1*(x2 - x1^2) - 2*(1 - x1),200*(x2 - x1^2))
}
m <- optim(c(-1.2,1),fr,control=list(trace=TRUE))
mg <- optim(c(-1.2,1),fr,gr=gr,method="BFGS",
  control=list(trace=TRUE))

# ------------------------------------------------------------------
## Nelinearna optimizacija / OECD

oecd <- read.table(file="OECD.dat",header=TRUE,sep="",row.names=1)
pairs(oecd)
attach(oecd)
plot(agr,pcinc,pch="+")

# linearna regresija
lin <- lm(pcinc ~ agr)
abline(lin,col="green")
lp <- lin$coef[2]*agr + lin$coef[1]
sum((lp - pcinc)^2)

# eksponentna z linearno regresijo
pi <- log(pcinc); m <- lm(pi ~ agr)
b <- exp(m$coef[1]); a <- exp(m$coef[2])
pl <- function(x){b*a^x}
points(agr,pl(agr),col="red",pch=16)

# metoda najmanjsih kvadratov
f <- function(t,p){a <- p[1]; b <- p[2]; b*a^t}
E <- function(p){d <- f(agr,p) - pcinc; sum(d^2)}
p0 <- c(a,b); best <- optim(p0,E)
E(p0)
best
pr <- function(x){f(x,best$par)}
points(agr,pr(agr),col="blue",pch=16)
d <- seq(0,84,2); lines(spline(d,pr(d)),col="blue")

# ------------------------------------------------------------------