## source("eigen.R")

## non positive semidefinite matrix
A <- cbind( c(1,2), c(2,3) )
print(A)
print(eigen(A))

## correlation matrix (2x2)
rho <- 0.6
A <- (1-rho)*diag(2) + rho
print(A)
print(eigen(A))

## correlation matrix (5x5)
rho <- 0.6
A <- (1-rho)*diag(5) + rho
print(A)
print(eigen(A))

library(ellipse)
## iso-distance curve, x'x=1
plot(ellipse(0),type="l",xlim=c(-5,5),ylim=c(-5,5))

## iso-distance curve, x'Ax=1, A correlation matrix with rho=0.8
plot(ellipse(0.8),type="l",xlim=c(-5,5),ylim=c(-5,5))

## iso-distance curve, x'Ax=1,
## A <- cbind( c(1,1.6), c(1.6,4))
plot(ellipse(0.8,scale=c(1,2)),type="l",xlim=c(-5,5),ylim=c(-5,5))

library(mvtnorm)

x1 <- x2 <- seq(-5, 5, length= 50)

nd <- function(x1,x2,rho) {
  A <- (1-rho)*diag(2) + rho
  r <- dmvnorm(cbind(x1,x2), mean=c(0,0), sigma=A)
  return(r)
}

## bivariate normal density, perspective and contour plots
rho <- -0.7
par(mfrow=c(1,2))
z <- outer(x1, x2, nd, rho=rho)
persp(x1, x2, z, theta=30, phi=30, expand=0.5, col="lightblue",
      ltheta=120, shade = 0.75)
contour(x1, x2, z,nlevels=10)
##points(rmvnorm(1000,sigma=(1-rho)*diag(2) + rho),pch=20,cex=.5,col="cyan")
par(mfrow=c(1,1))