I found an easy way to do it using trans3d() and the transformation matrix.
Maybe in the future somebody would ask about it:

# For plotting the bivariate gaussian density
library(mvtnorm)

x<-seq(-5,5,by=0.1)
y<-seq(-3,7,by=0.1)

# Joint density
f<-function(x,y)
dmvnorm(cbind(x,y),mean=c(0,2),sigma=matrix(c(3,2,2,4),ncol=2))

# Marginal densities
fx<-function(x) dnorm(x,0,3)
fy<-function(y) dnorm(y,2,4)

z<-outer(x,y,f)

# Store the transformation matrix
mat<-persp(x,y,z,theta=30,phi=30,zlim=c(0,0.15))

# Draw the marginals
l1<-trans3d(rep(min(x),length(x)),y,fy(y),pmat=mat)
l2<-trans3d(x,rep(max(y),length(y)),fx(x),pmat=mat)
lines(l1$x,l1$y,lwd=1)
lines(l2$x,l2$y,lwd=1)

	[[alternative HTML version deleted]]

