## Code for Stat 624 lecure on constrained optimization ## Written by Robert Richardson # Sample from mixture of normals set.seed(12345) dat <- rnorm(100000,sample(c(-10,0,10),1000,replace=TRUE,prob=c(1,2,3)),sample(c(4.5,1.75,.6),100,replace=TRUE,prob=c(1,2,3))) hist(dat,breaks = 30) log_lik <- function(pars){ mus <- pars[1:3] sigs <- pars[4:6] weights <- pars[7:9] -sum(log(weights[1]*dnorm(dat,mus[1],sigs[1])+ weights[2]*dnorm(dat,mus[2],sigs[2])+ weights[3]*dnorm(dat,mus[3],sigs[3]))) } mu <- c(-5.0,0.0,5.0) sigs <- c(3.0,3.0,3.0) weights <- c(1.0,1.0,1.0)/3 opt <- optim(c(mu,sigs,weights),log_lik) ### Impose constaints in the likelihood log_lik <- function(pars){ mus <- pars[1:3] sigs <- pars[4:6] weights <- c(pars[7:8],1-sum(pars[7:8])) if(any(weights < 0) | sum(weights > 1)){ return(100000000000) } else { -sum(log(weights[1]*dnorm(dat,mus[1],sigs[1])+ weights[2]*dnorm(dat,mus[2],sigs[2])+ weights[3]*dnorm(dat,mus[3],sigs[3]))) } } mu <- c(-5.0,0.0,5.0) sigs <- c(1.0,1.0,1.0) weights <- c(.1,.1) opt <- optim(c(mu,sigs,weights),log_lik) mus <- opt$par[1:3] sigs <- opt$par[4:6] weights <- c(opt$par[7:8],1-sum(opt1$par[7:8])) hist(dat,freq=FALSE,breaks=30) f <- function(x){ weights[1]*dnorm(x,mus[1],sigs[1]) + weights[2]*dnorm(x,mus[2],sigs[2])+ weights[3]*dnorm(x,mus[3],sigs[3]) } curve(f,from=-40,to=40,add=TRUE,col='red') ## Broyden-Fletcher-Goldfarb-Shanno!! b.lower <- c(-Inf,-Inf,-Inf,0,0,0,0,0) b.upper <- c(Inf,Inf,Inf,Inf,Inf,Inf,1,1) inits <- c(-10.0,0.0,10.0,1.0,1.0,1.0,0.3,0.3) opt1 <- optim(inits,log_lik,method="L-BFGS-B",lower=b.lower,upper=b.upper,control=list(maxit=1000)) opt1 mus <- opt1$par[1:3] sigs <- opt1$par[4:6] weights <- c(opt1$par[7:8],1-sum(opt1$par[7:8])) hist(dat,freq=FALSE,breaks=30) f <- function(x){ weights[1]*dnorm(x,mus[1],sigs[1]) + weights[2]*dnorm(x,mus[2],sigs[2])+ weights[3]*dnorm(x,mus[3],sigs[3]) } curve(f,from=-40,to=40,add=TRUE,col='red') ## Linear Programming log_lik <- function(pars){ mus <- pars[1:3] sigs <- pars[4:6] weights <- c(pars[7:8],1-sum(pars[7:8])) -sum(log(weights[1]*dnorm(dat,mus[1],sigs[1])+ weights[2]*dnorm(dat,mus[2],sigs[2])+ weights[3]*dnorm(dat,mus[3],sigs[3]))) } inits <- c(-10.0,0.0,10.0,1.0,1.0,1.0,0.3,0.3) # Optimize subject to constraints ui %*% pars - ci >= 0 ui <- rbind(cbind(matrix(0,5,3),diag(5)), cbind(matrix(0,2,6),-diag(2)), c(rep(0,6),-1,-1)) ci <- c(rep(0,5),-1,-1,-1) ui%*%inits - ci opt2 <- constrOptim(inits,log_lik,grad=NULL,ui=ui,ci=ci) mus <- opt2$par[1:3] sigs <- opt2$par[4:6] weights <- c(opt2$par[7:8],1-sum(opt1$par[7:8])) hist(dat,freq=FALSE,breaks=30) f <- function(x){ weights[1]*dnorm(x,mus[1],sigs[1]) + weights[2]*dnorm(x,mus[2],sigs[2])+ weights[3]*dnorm(x,mus[3],sigs[3]) } curve(f,from=-40,to=40,add=TRUE,col='red')