## Perspective plot -- unbiased estimate of variance
res <- 37
x <- y <- seq(0, 1, len=res)
Z <- matrix(NA, nrow=res, ncol=res)
calcW <- function(x, y) {
  if (y==1) return(c(0, 0, 1))
  w3 <- y
  d <- y/2
  w2 <- (x-d)/(1-2*d)*(1-w3)
  w1 <- 1 - w2 - w3
  return(c(w1, w2, w3))
}
for (i in 1:res) {
  for (j in 1:res) {
    w <- calcW(x[i], y[j])
    if (any(w < 0 | w > 1)) Z[i,j] <- 0
    else Z[i,j] <- cov.wt(matrix(c(1,3,2),ncol=1), wt=w)$cov
  }
}
pobj <- persp(x, y, Z, box=FALSE, theta=150, axes=FALSE, lwd=0.5, phi=25, d=5, border="blue")
points(trans3d(0.5,1,0.52,pmat=pobj),pch=19)
text(trans3d(0.5,1,0.52,pmat=pobj)$x+.03,trans3d(0.5,1,0.52,pmat=pobj)$y,labels="2")
points(trans3d(0,0,0.52,pmat=pobj),pch=19)
text(trans3d(0,0,0.52,pmat=pobj)$x+.03,trans3d(0,0,0.52,pmat=pobj)$y,labels="3")
points(trans3d(1,0,0.52,pmat=pobj),pch=19)
text(trans3d(1,0,0.52,pmat=pobj)$x-.03,trans3d(1,0,0.52,pmat=pobj)$y,labels="1")
points(trans3d(0.5,1/3,var(1:3),pmat=pobj),pch=19)
text(trans3d(0.5, 1/3, var(1:3), pmat=pobj)$x-.03, trans3d(0.5, 1/3, var(1:3), pmat=pobj)$y,labels=expression(hat(w)))
x1 <- x; y1 <- y; Z1 <- Z

## Perspective plot -- plugin variance
res <- 37
x <- y <- seq(0, 1, len=res)
Z <- matrix(NA, nrow=res, ncol=res)
calcW <- function(x, y) {
  if (y==1) return(c(0, 0, 1))
  w3 <- y
  d <- y/2
  w2 <- (x-d)/(1-2*d)*(1-w3)
  w1 <- 1 - w2 - w3
  return(c(w1, w2, w3))
}
for (i in 1:res) {
  for (j in 1:res) {
    w <- calcW(x[i], y[j])
    if (any(w < 0 | w > 1)) Z[i,j] <- 0
    else Z[i,j] <- cov.wt(matrix(c(1,3,2),ncol=1), wt=w, method="ML")$cov
  }
}
pobj <- persp(x, y, Z, box=FALSE, theta=150, axes=FALSE, lwd=0.5, phi=25, d=5, border="blue")
points(trans3d(0.5,1,0,pmat=pobj),pch=19)
text(trans3d(0.5,1,0,pmat=pobj)$x+.03,trans3d(0.5,1,0,pmat=pobj)$y,labels="2")
points(trans3d(0,0,0,pmat=pobj),pch=19)
text(trans3d(0,0,0,pmat=pobj)$x+.06,trans3d(0,0,0,pmat=pobj)$y,labels="3")
points(trans3d(1,0,0,pmat=pobj),pch=19)
text(trans3d(1,0,0,pmat=pobj)$x-.03,trans3d(1,0,0,pmat=pobj)$y,labels="1")
v <- cov.wt(matrix(c(1,3,2),ncol=1), method="ML")$cov
points(trans3d(0.5,1/3,v,pmat=pobj),pch=19)
text(trans3d(0.5, 1/3, v, pmat=pobj)$x-.03, trans3d(0.5, 1/3, v, pmat=pobj)$y,labels=expression(hat(w)))
x2 <- x; y2 <- y; Z2 <- Z

## RGL
require(rgl)
surface3d(x1, y1, Z1, col="gray80")
surface3d(x2, y2, Z2, col="gray80")