Efficent way of creating a contour

48 views
Skip to first unread message

Brian Weaver

unread,
Oct 5, 2010, 3:23:20 PM10/5/10
to ggp...@googlegroups.com
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)) )


constrained_optim.txt
contour_failuretime_quantile.pdf

Brandon Hurr

unread,
Oct 5, 2010, 3:57:37 PM10/5/10
to Brian Weaver, ggp...@googlegroups.com
I had a discussion with Baptiste a while back about implementing patternGrob in the gridExtra package as a fill for geom_rect(), but it was well over my head. Perhaps you could use it to fill your curve area instead of drawing the lines by hand? It will probably take you many more hours to figure out how though... 

I'd be willing to help if I can, but I looked into it and just implementing patternGrob as a fill for geom_bar() was way over my head. 

If it's any consolation, that is a very nice graph. Even if I don't know what it means. 

Brandon



--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: http://gist.github.com/270442
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

Hadley Wickham

unread,
Oct 6, 2010, 10:13:53 PM10/6/10
to Brian Weaver, ggp...@googlegroups.com
Couldn't you automate the production of the lines? Use contourLines
to calculate the contours, then pick the appropriate level, and use
some basic geometry to work out where evenly spaced lines would
intersect with the curve. It would probably take a couple of hours to
work out the details, but once you had it, it would save a lot of time
for future plots.

Hadley

> --
> You received this message because you are subscribed to the ggplot2 mailing
> list.
> Please provide a reproducible example: http://gist.github.com/270442
>
> To post: email ggp...@googlegroups.com
> To unsubscribe: email ggplot2+u...@googlegroups.com
> More options: http://groups.google.com/group/ggplot2
>

--
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

Brian Weaver

unread,
Oct 7, 2010, 10:15:15 AM10/7/10
to ggp...@googlegroups.com
I had considered that but the issue is still having the cross hatch lines stop at the contours, which I would imagine would still take a lot of time picking those points, especially as the situation changes. What I was envisioning was first creating the plot that contained nothing but the cross hatched lines (which could be automated).  Then putting on top of these lines the standard error contours that had the color white in between the bins and then finally putting the cost contours on last.  In general what I am seeing is creating the plot layer by layer.  I have not had the time to figure this out yet, however, but if I ever get that time I will resend with the code.

Thank you all for your time and help.

Brian
Reply all
Reply to author
Forward
0 new messages