set.seed(3223)
Sigma<-rWishart(1, 3, diag(c(1,1)))
Sigma<-Sigma[,,1]

Prec<-solve(Sigma)

library(MASS)
X<-mvrnorm(1000,c(0,0),Sigma)

plot(X)
plot(X,pch=19)
plot(X,pch=15)
plot(X,pch=3)
plot(X,pch=3,cex=0.5)

help(points)

xv<-yv<-seq(-4.5,4.5,by=0.1)

MVNdens<-function(xval,yval,Pmat){

	tval<--0.5*(Pmat[1,1]*xval^2+2*Pmat[1,2]*xval*yval+Pmat[2,2]*yval^2)
	
	tval<--2*0.5*log(2*pi)*0.5*log(det(Pmat))+tval
	
	return(tval)

}

MVNdens(1,2,Prec)

Dmat<-outer(xv,yv,MVNdens,Pmat=Prec)

image(xv,yv,t(Dmat))
contour(xv,yv,t(Dmat),add=T)
points(X,pch=19,cex=0.5)

image(xv,yv,t(Dmat))
contour(xv,yv,t(Dmat),add=T)

pdf("NormPlot.pdf")
image(xv,yv,t(Dmat))
contour(xv,yv,t(Dmat),add=T)
dev.off()

range(Dmat)
pdf("NormPlot2.pdf")
image(xv,yv,t(Dmat))
contour(xv,yv,t(Dmat),add=T,levels=seq(-4,1,by=0.5))
dev.off()

pdf("NormPlot2.pdf",paper="USr",width=11,height=9)
par(pty="s")
image(xv,yv,t(Dmat))
contour(xv,yv,t(Dmat),add=T,levels=seq(-4,1,by=0.5))
dev.off()


help(pdf)