options(OutDec = ",") 
#==============================================================================
# Carregando funcoes
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(mvtnorm))

#==============================================================================
# Analisando graficos da normal bivariada
#==============================================================================
nx  <- 100
ny  <- 100
x   <- seq(-4,4,length=nx)
y   <- seq(-4,4,length=ny)
z   <- matrix(NA,nx,ny)
mu  <- c(0,0)
S   <- matrix(c(1, .7, .7, 1),2,2)
for(i in 1:nx)
    for(j in 1:ny)
        z[i,j] <- dmvnorm(c(x[i],y[j]), mean=mu,sigma=S)

#==============================================================================
# Grafico 1: funcao de densidade de probabilidade
#==============================================================================
persp(x,y,z, phi = 45, theta = 30,col="magenta",xlab=expression(x),
      ylab=expression(y),zlab=expression(f(x,y)))

#==============================================================================
# Grafico 2: curvas de contorno da funcao de densidade de probabilidade
#==============================================================================
contour(x,y,z,col="red",lwd=2,xlab=expression(x),ylab=expression(y),bty="n")

#==============================================================================
# Gerando dados da normal bivariada
# Var(X) = 1; Var(Y) = 1 e Cov(X,Y) = 0,7
#==============================================================================
set.seed(5489) # fixa a sequencia de numeros gerados
print(mu)      # Vetor de medias
print("====================")
print(S)       # Matriz de covariancias
print("====================")
w   <- rmvnorm(100000,mean=mu,sigma=S)
print(head(w))
print("====================")

[1] 0 0
[1] "===================="
     [,1] [,2]
[1,]  1,0  0,7
[2,]  0,7  1,0
[1] "===================="
             [,1]        [,2]
[1,] -0,867344397 -0,75834358
[2,] -0,090669610  0,56366994
[3,]  2,002714275  1,37908761
[4,] -0,799916903 -0,25788300
[5,] -0,068928331  0,51429148
[6,] -0,004930412  0,02791009
[1] "===================="

#==============================================================================
# Estatisticas descritivas
#==============================================================================
apply(w,2,"mean")  # Media amostral
print(cov(w))      # Covariancia amostral
print("====================")
print(cor(w))      # Correlacao amostral
print("====================")

          [,1]      [,2]
[1,] 0,9968351 0,6982036
[2,] 0,6982036 0,9967862
[1] "===================="
          [,1]      [,2]
[1,] 1,0000000 0,7004375
[2,] 0,7004375 1,0000000
[1] "===================="

#==============================================================================
# Grafico 3: dispersao
#==============================================================================
plot(w,pch=".",main="",xlab=expression(x),ylab=expression(y),bty="n")
contour(x,y,z,col="blue",lwd=2,lty=2,xlab=expression(x),ylab=expression(y),
    add=T)

#==============================================================================
# Grafico 4: histograma de x
#==============================================================================
par(bty="n",lab=c(10,10,5))
hist(w[,1],prob=T,nclass=50,main="",xlab=expression(x),ylab="Densidade",
    xlim=c(-5,5),ylim=c(0,0.45))
lines(x,dnorm(x),lwd=2,col="red")

#==============================================================================
# Grafico 4: histograma de x
#==============================================================================
par(bty="n",lab=c(10,10,5))
hist(w[,2],prob=T,nclass=50,main="",xlab=expression(x),ylab="Densidade",
    xlim=c(-5,5),ylim=c(0,0.45))
lines(y,dnorm(y),lwd=2,col="red")