All,
I have been working on a constrained optimization problem where I find the solution graphically (think Lagrange multipliers). The idea is I want to find the number of items (n) and the number of measurements per item (m) that minimizes some cost function subject to the constraint that the estimation precision of a function of interest (a standard error) is less than a certain value. So the idea behind the plot is to have a contour that shows the contours of the cost function, the contours of the standard error, and a region that is not allowed since these values of (n,m) do not satisfy the constraint. My professor wants me to "cross hatch out" the region that is not allowed. Attached to this email is one of these plots so that you can better understand my description. The dashed lines represent the contours for the standard error and the thick solid contours represent the cost function. Also the asterisk is the optimum point. I can do this plot by hand, but, one plot takes around 1-2 hours to create. My question is is anyone familiar of a quicker way of doing this? Currently the way I do it is create the contours and then add a geom_line for each individual line in the cross hatch region. My code is below and I am attaching a data set that was used to create the plot. I thank you all for your help ahead of time.
Brian
library(ggplot2)
dat <- read.table("your directory here...constrained_optim.txt",header = T)
#Selecting the values that satisfy the constraint and saving as a separate object
t.id <- which(dat$SE <= 0.8)
xx <- dat[
t.id,]
tt <- which(xx[,3] == min(xx[,3]))
xx[tt,] ##Solution to optimum problem
##Creating the plot object
v <- ggplot(dat,aes(x = n, y = m, z = cost))
p <- v + geom_contour(breaks = seq(48600,80000,by = 7500),size = 1 ,color = "black" )+
stat_contour(data = dat,aes(x = n,y = m, z = SE),
lty = 2, breaks = c(0.7906603,0.75,0.7)) +
geom_point(aes(x = dat$n, y = dat$m), size =0) +
geom_point(aes(x = xx[tt,1],y = xx[tt,2]),shape = 8,size = 4)+theme_bw() +
ylab("Number of Measurements (m)") + xlab("Number of Units (n)")+
geom_text(aes(x = 11.3, y = 20),label = 49,angle = 45, size = 4, hjust = 0)+
geom_text(aes(x = 13.8, y = 20),label = 56,angle = 45, size = 4, hjust = 0)+
geom_text(aes(x = 16.3, y = 20),label = 64,angle = 45, size = 4, hjust = 0)+
geom_text(aes(x = 18.7, y = 20),label = 71,angle = 45, size = 4, hjust = 0)+
geom_text(aes(x = 20, y = 4.4),label = 0.79, size = 3, hjust = 0)+
geom_text(aes(x = 20, y = 5.6),label = 0.75, size = 3, hjust = 0)+
geom_text(aes(x = 20, y = 8.8),label = 0.7, size =3, hjust = 0)
p
#Now adding the cross hatch region
pp <- p+geom_line(aes(x = c(19,20),y = c(3,4)) )+
geom_line(aes(x = c(18,19.5),y = c(3,4.6)) )+
geom_line(aes(x = c(17,18.9),y = c(3,5)) )+
geom_line(aes(x = c(16,18.4),y = c(3,5.4)) )+
geom_line(aes(x = c(15,17.8),y = c(3,5.8)) )+
geom_line(aes(x = c(14,17.3),y = c(3,6.4)) )+
geom_line(aes(x = c(13,16.8),y = c(3,6.9)) )+
geom_line(aes(x = c(12,16.4),y = c(3,7.4)) )+
geom_line(aes(x = c(11,16),y = c(3,8)) )+
geom_line(aes(x = c(11,15.7),y = c(4,8.7)) )+
geom_line(aes(x = c(11,15.4),y = c(5,9.5)) )+
geom_line(aes(x = c(11,15),y = c(6,10)) )+
geom_line(aes(x = c(11,14.8),y = c(7,10.8)) )+
geom_line(aes(x = c(11,14.6),y = c(8,11.6)) )+
geom_line(aes(x = c(11,14.4),y = c(9,12.4)) )+
geom_line(aes(x = c(11,14.2),y = c(10,13.2)) )+
geom_line(aes(x = c(11,14),y = c(11,14)) )+
geom_line(aes(x = c(11,13.9),y = c(12,14.9)) )+
geom_line(aes(x = c(11,13.7),y = c(13,15.7)) )+
geom_line(aes(x = c(11,13.6),y = c(14,16.5)) )+
geom_line(aes(x = c(11,13.5),y = c(15,17.5)) )+
geom_line(aes(x = c(11,13.3),y = c(16,18.4)) )+
geom_line(aes(x = c(11,13.3),y = c(17,19.3)) )+
geom_line(aes(x = c(11,13),y = c(18,20)) )+
geom_line(aes(x = c(11,12),y = c(19,20)) )
pp
##Positive sloped lines end here
pp+geom_line(aes(x = c(13,13.2),y = c(20,19.8)) )+
geom_line(aes(x = c(12,13.3),y = c(20,18.6)) )+
geom_line(aes(x = c(11,13.5),y = c(20,17.5)) )+
geom_line(aes(x = c(11,13.6),y = c(19,16.3)) )+
geom_line(aes(x = c(11,13.8),y = c(18,15.1)) )+
geom_line(aes(x = c(11,14),y = c(17,14)) )+
geom_line(aes(x = c(11,14.3),y = c(16,12.8)) )+
geom_line(aes(x = c(11,14.6),y = c(15,11.5)) )+
geom_line(aes(x = c(11,15.1),y = c(14,9.9)) )+
geom_line(aes(x = c(11,16),y = c(13,8)) )+
geom_line(aes(x = c(11,20),y = c(12,3)) )+
geom_line(aes(x = c(11,19),y = c(11,3)) )+
geom_line(aes(x = c(11,18),y = c(10,3)) )+
geom_line(aes(x = c(11,17),y = c(9,3)) )+
geom_line(aes(x = c(11,16),y = c(8,3)) )+
geom_line(aes(x = c(11,15),y = c(7,3)) )+
geom_line(aes(x = c(11,14),y = c(6,3)) )+
geom_line(aes(x = c(11,13),y = c(5,3)) )+
geom_line(aes(x = c(11,12),y = c(4,3)) )