# Optimization
# October 22, 2015
# Sunthud Pornprasertmanit


#############################################
# Single-variable problems
############################################

### Example 1

FUN1 <- function(x) x * log(x) + x^2 - 5

uniroot(FUN1, c(0.00001, 100))

plot(seq(0.01, 5, 0.01), 
   FUN1(seq(0.01, 5, 0.01)), type = "l")
abline(h = 0, lty = 2)

FUN1e <- function(x) x^3 - 3*x^2 - 18*x + 40
plot(seq(-5, 7, 0.01), 
   FUN1e(seq(-5, 7, 0.01)), type = "l")
uniroot(FUN1e, c(-5, 7))
   
###########################################

### Example 2

FUN2 <- function(intcept, pmiss) {
	logx <- intcept + 0.7*0.5 
	1/(1 + exp(-logx)) - pmiss
}
uniroot(FUN2, c(-5, 5), pmiss = 0.2)

plot(seq(-5, 5, 0.01), 
   FUN2(seq(-5, 5, 0.01), pmiss = 0.2) + 0.2, type = "l")

#############################################

### Example 3

FUN3 <- function(n, pow, d, alpha = 0.05) {
	crit <- qt(1 - alpha/2, df = n - 2)
	obspow <- 1 - pt(crit, df = n - 2, 
		ncp = d * sqrt(n/4))
	obspow - pow
}
uniroot(FUN3, c(5, 1000), d = 0.5, pow = 0.8)
plot(seq(5, 200, 1), 
   FUN3(seq(5, 200, 1), d = 0.5, pow = 0.8) + 0.8, type = "l")
abline(h = 0.8, lty = 2)

############################################

### Example 4

x <- attitude$rating

sum(dchisq(x, df = 50, log = TRUE))
sum(dchisq(x, df = 60, log = TRUE))
sum(dchisq(x, df = 70, log = TRUE))
sum(dchisq(x, df = 80, log = TRUE))
sum(dchisq(x, df = 90, log = TRUE))

x <- rchisq(200, 10)

FUN4 <- function(x, dat) {
	(2*1) -2 * sum(dchisq(dat, df = x, log = TRUE))
}
nlminb(5, FUN4, dat = x, lower = 2, upper = 1000)
optim(5, FUN4, dat = x, lower = 2, upper = 1000, method = "L-BFGS-B", hessian = TRUE)
fit <- optim(5, FUN4, dat = x, lower = 2, upper = 1000, method = "L-BFGS-B", hessian = TRUE)
infomat <-solve(fit$hessian)
se <- sqrt(diag(fisher_info))
crit <- qnorm(0.975)
c(fit$par - crit*se, fit$par + crit*se)

# Local vs Global minima

FUN41 <- function(x) x^4 + x^3 - 3*x^2 + x + 2

plot(seq(-3, 3, 0.1), FUN41(seq(-3, 3, 0.1)), type = "l")

nlminb(2, FUN41)
nlminb(-2, FUN41)


###################################

### Example 5

findload <- function(k, target) {
	loading <- k * seq(0.5, 0.9, 0.1)
	errvar <- 1 - (loading^2)
	omega <- sum(loading)^2/(sum(loading)^2 + sum(errvar))
	omega - target
}

uniroot(findload, c(0.1, 1.5), target = 0.8)

#####################################
# Multi-variable problems
#####################################

### Example 1

MUL1 <- function(vec) {
	rhs1 <- vec[1]^vec[2]
	rhs2 <- vec[1] * vec[2]
	(rhs1 - 5)^2 + (rhs2 - 100)^2
}
nlminb(c(5, 2), MUL1)

MUL1e <- function(vec) {
	rhs1 <- log(vec[1] + vec[2])
	rhs2 <- vec[1]/vec[2]
	(rhs1 - (-100))^2 + (rhs2 - 5)^2
}
nlminb(c(9, 3), MUL1e)

MUL1s <- function(vec) {
	rhs1 <- vec[1]*vec[2]
	rhs2 <- vec[1]/vec[2]
	(rhs1 - 20)^2 + (rhs2 - 5)^2
}
nlminb(c(9, 3), MUL1s)
nlminb(c(-9, -3), MUL1s)

################################################

### Example 2

MUL2 <- function(par, vals, pm) {
	diff1 <- log(pm[1]/(1 - pm[1])) - par[1] - par[2]*vals[1] 
	diff2 <- log(pm[2]/(1 - pm[2])) - par[1] - par[2]*vals[2]
	diff1^2 + diff2^2
}
nlminb(c(0, 0), MUL2, vals = c(0, 1), pm = c(0.4, 0.2))

##############################################

### Example 3

findn1 <- function(n2, pow = 0.8, d = 0.5, alpha = 0.05) {
	FUN <- function(n1, pow, d, alpha, n2) {
		crit <- qt(1 - alpha/2, df = n1 + n2 - 2)
		obspow <- 1 - pt(crit, df = n1 + n2 - 2, 
			ncp = d * sqrt(n1*n2/(n1 + n2)))
		obspow - pow
	}
	uniroot(FUN, c(5, 1000), pow = pow, d = d, alpha = alpha, n2 = n2)$root
}

plot(apply(matrix(35:250), 1, findn1), 35:250, type = "l", xlab = "n1", ylab = "n2", main = "power = .80")


	
MUL3 <- function(ns, pow, d, alpha = 0.05) {
	crit <- qt(1 - alpha/2, df = ns[1] + ns[2] - 2)
	obspow <- 1 - pt(crit, df = ns[1] + ns[2] - 2, 
		ncp = d * sqrt(ns[1]*ns[2]/(ns[1] + ns[2])))	
	cost <- (2*ns[1] + ns[2])/100000
	(obspow - pow)^2 + cost
}

out <- nlminb(c(20, 10), MUL3, d = 0.5, pow = 0.8, lower = c(2,2))
ns <- out$par
crit <- qt(1 - 0.05/2, df = ns[1] + ns[2] - 2)
obspow <- 1 - pt(crit, df = ns[1] + ns[2] - 2, 
	ncp = 0.5 * sqrt(ns[1]*ns[2]/(ns[1] + ns[2])))
cost <- (2*ns[1] + ns[2])
c(ns, obspow, cost)


MUL32 <- function(n2, pow, d, alpha = 0.05) {
	FUN <- function(n1, pow, d, alpha, n2) {
		crit <- qt(1 - alpha/2, df = n1 + n2 - 2)
		obspow <- 1 - pt(crit, df = n1 + n2 - 2, 
			ncp = d * sqrt(n1*n2/(n1 + n2)))
		obspow - pow
	}
	out <- uniroot(FUN, c(5, 1000), pow = pow, d = d, alpha = alpha, n2 = n2)
	n1 <- out$root
	(2*n1 + n2)
}
nlminb(50, MUL32, d = 0.5, pow = 0.8, lower = 2)

#####################################

### Example 4

x <- attitude$rating

dnorm(x, mean = 0, sd = 1, log = TRUE)
sum(dnorm(x, mean = 0, sd = 1, log = TRUE))

MUL4 <- function(x, dat) {
	-sum(dnorm(dat, mean = x[1], sd = x[2], log = TRUE))
}
nlminb(c(0, 1), MUL4, dat = x, lower = c(-Inf, 0.001))

nlminb(c(1000, 200), MUL4, dat = x, lower = c(-Inf, 0.001))
nlminb(c(-20, 5), MUL4, dat = x, lower = c(-Inf, 0.001))
nlminb(c(2500, 0.001), MUL4, dat = x, lower = c(-Inf, 0.001))
nlminb(c(-10000, 1000), MUL4, dat = x, lower = c(-Inf, 0.001))
nlminb(c(50, 10), MUL4, dat = x, lower = c(-Inf, 0.001))

#############################################

### Example 5

MUL5 <- function(x, target) {
	loading <- x
	errvar <- 1 - (loading^2)
	omega <- sum(loading)^2/(sum(loading)^2 + sum(errvar))
	(omega - target)^2
}
nlminb(c(0.9, 0.45, 0.7, 0.5, 0.5), MUL5, target = 0.8, lower = 0.4, upper = 0.95)
nlminb(c(0.45, 0.45, 0.45, 0.5, 0.9), MUL5, target = 0.8, lower = 0.4, upper = 0.95)

#################################################

### Example 6

library(lavaan)
dat <- attitude[2:7]
colnames(dat) <- c("y1", "y2", "y3", "y4", "y5", "y6")
MUL6 <- function(w, dat) {
	w <- c(w, 1 - sum(w))
	p <- ncol(dat)
	script <- "f1 =~ y1 + y2 + y3 + y4 + y5 + y6"
	out <- lavaan::cfa(script, data = dat, std.lv = TRUE)
	LY <- inspect(out, "coef")$lambda
	TH <- inspect(out, "coef")$theta
	IMPLIED <- LY %*% t(LY)
	num <- t(w) %*% (IMPLIED) %*% w
	denom <- t(w) %*% (IMPLIED + TH) %*% w
	1 - num/denom
}

nlminb(rep(1, 5), MUL6, dat = dat)


#########################################
### Extra: Solving Simple Regression Problem
#########################################

x <- rnorm(100, 0, 1)
y <- 0.7*x + rnorm(100, 0, 0.5)

summary(lm(y ~ x))

FUNS <- function(param, x, y) {
	err <- y - (param[1] + param[2]*x)
	- sum(dnorm(err, mean = 0, sd = param[3], log = TRUE))
}
fit <- optim(c(0, 0.7, 0.5), FUNS, x = x, y = y, lower = c(-Inf, -Inf, 0.0001), method = "L-BFGS-B", hessian = TRUE)

fisher_info<-solve(fit$hessian)
prop_sigma<-sqrt(diag(fisher_info))