par(mfrow=c(2,2))

M <- 2
ngrd <- 50
l <- (2*ngrd+1)
lplt <- (M+1)*l^M

plt  <- scan("npb.big_down.plt")
mu   <- plt[seq(from=1,to=M)]
sig  <- plt[seq(from=M+1,to=M+M*M)]
zinc <- plt[seq(from=M+M*M+1,to=M+M*M+M)]
plt  <- plt[seq(from=M+M*M+M+1,to=M+M*M+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
x <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
y <- plt[seq(from=2,to=lplt,by=l*(M+1))]

persp(x,y,h,xlab="y_1",ylab="y_2",zlab="")
title("(a)")

Eh <- -sort(-cumsum(sort(zinc[1]*zinc[2]*h)))
Oh <- order(-h)
qh <- c(h[Oh[length(Eh[Eh>.95])]], h[Oh[length(Eh[Eh>.90])]],
        h[Oh[length(Eh[Eh>.75])]], h[Oh[length(Eh[Eh>.50])]],
        h[Oh[length(Eh[Eh>.25])]], h[Oh[length(Eh[Eh>.10])]])

contour(x,y,h,labex=0,xlab="y_1",ylab="y_2",levels=qh)
title("(b)")

one <- seq(1,1,length=l)

hx <- zinc[2]*h%*%one
hy <- zinc[1]*one%*%h
nx <- dnorm(x,mean=mu[1],sd=sqrt(sig[1]))
ny <- dnorm(y,mean=mu[2],sd=sqrt(sig[3]))

plot(c(y,y),c(hy,ny),type='n',xlab="y_2",ylab="")
lines(y,hy,lty=1)
lines(y,ny,lty=2)
title("(c)")

plot(c(x,x),c(hx,nx),type='n',xlab="y_1",ylab="")
lines(x,hx,lty=1)
lines(x,nx,lty=2)
title("(d)")

