[R] Draw level lines on the surface of a bivariate function
Duncan Murdoch
murdoch at stats.uwo.ca
Sat Mar 11 19:03:06 CET 2006
On 3/11/2006 12:17 PM, Cuvelier Etienne wrote:
> Hello,
> Is it possible to draw level lines on the surface of a bivariate function?
>
> In the following example, to draw surface and levels lines for a
> multivariate normal law,
> I use persp, trans3d, contourLines and lines,
> but if the lines are correctly drawn, some parts of them are, of
> course, visible
> even if they are drawn on a non visible "face".
>
> Any suggestion to avoid this problem ?
You could try to work out which parts of the contour should be hidden,
based on their direction---but that's probably hard.
You could try the rgl package, but you'll need some fiddling, because
when you draw lines at the same place you've drawn a surface, it's
mostly rounding error that determines which one shows up. rgl also
doesn't do the automatic rescaling that persp does. However, this code
does something pretty close to what you want:
library(rgl)
open3d()
bg3d("white")
surface3d(x,y,20*z,col='red')
for (i in 1:length(lev))
lines3d(lev[[i]]$x, lev[[i]]$y, 20*lev[[i]]$level+0.01, col='black')
bbox3d(back='lines',xat=pretty(range(x)), yat=pretty(range(y)),
zat=pretty(range(z))*20, zlab=pretty(range(z)))
Duncan Murdoch
> Thank you
>
> Etienne
>
> Example :
>
> trans3d <- function(x,y,z, pmat)
> {
> tr <- cbind(x,y,z,1) %*% pmat
> list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4])
> }
>
> mean1 = c(2,4)
> cov1 = matrix(c(0.7, 0.2,0.2,0.7), ncol=2)
>
> x = seq(0, 5, by = 0.1)
> y = seq(0,10, by= 0.1)
> z = matrix(nrow = length(x), ncol=length(y))
>
> for(i in 1:length(x))
> for(j in 1:length(y))
> {
> z[i,j]=dmvnorm(c(x[i],y[j]),mean1,cov1)
> }
>
> pmat1 = persp(x, y, z, col= "red",shade = 0.25, border=NA)
>
> lev =contourLines(x,y,z, nlevel = 10)
> for(i in 1:length(lev))
> lines(trans3d(lev[[i]]$x, lev[[i]]$y, lev[[i]]$level, pmat1), col="black")
>
>
> --
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list