Построение параметрического распределения многомерной нормали в R

Question

Учитывая следующую выборку многомерного нормального:

mu=rep(0,2)
Sigma=matrix(c(1,0,0,1),2,2)
require(MASS)
X=mvrnorm(n=100,mu,Sigma)

Я хотел бы вычислить это параметрическая функция плотности. Раньше это делалось так:

require(rgl)
require(mnormt)
zX=dmnorm(as.matrix(X),mean=colMeans(X),varcov=cov(X))
plot3d(X[,1],X[,2],zX) 

но это всплывает график точек (из RStudio):

Интересно, как я могу получить графические окна RStudio, которые также рисуют поверхность с этими точками? Это я мог бы сделать это с непараметрической плотностью:

ngrid = 50
Zbiksm = kde2d(X[,1],X[,2],n=ngrid)
Zbiksm_den = Zbiksm$z # valor de la densidad en cada punto del grid
persp(x=Zbiksm$x,y=Zbiksm$y,z=Zbiksm_den,xlab="x",ylab="y",
      zlab=expression(K(x,y)),theta=-35,axes=TRUE,box=TRUE)

но не знаю, как это сделать в этом случае.

В качестве альтернативы, это также помогает мне собрать изображения Surface3D и ScatterPlot.