In [1]:
# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Descricao
#
# Formas de gerar numeros pseudos-aleatorios da normal multivariada
#
#==============================================================================
set.seed(5485)                  # fixando a semente pseudo-aleatorio
p     <- 3                      # dimensao
n     <- 100000                 # tamanho da amostra
In [2]:
#==============================================================================
# estrutura de covariancia (correlacao)
#==============================================================================
SIGMA <- matrix(c(1,0.9,0.8,0.9,1,0.7,0.8,0.7,1),p,p,T)
print(SIGMA)
     [,1] [,2] [,3]
[1,]  1,0  0,9  0,8
[2,]  0,9  1,0  0,7
[3,]  0,8  0,7  1,0
In [3]:
#==============================================================================
# normais independentes
#==============================================================================
z     <- matrix(rnorm(p*n),n,p) # note a transposicao
In [4]:
#==============================================================================
# Utilizando a decomposicao SVD
#==============================================================================
S     <- svd(SIGMA)         # decomposicao de SIGMA
V     <- S$u                # matriz V de autovetores
Dhalf <- diag(sqrt(S$d))    # matriz D^{1/2} (raizes quadradas dos autovalores)
H     <- V%*%Dhalf          # construcao da matriz H
x     <- z%*%t(H)           # gerando valores (note a transposicao)

print("==== Médias ====")
print(apply(x,2,"mean"))  # vetor (linhas) de medias amostral

print("==== Covariâncias ====")
print(var(x))             # matriz de covariancias amostral

print("==== Correlações ====")
print(cor(x))             # matriz de correlacoes amostral

print("==== Covariâncias ====")
print(SIGMA)
[1] "==== Médias ===="
[1]  3,225724e-06 -3,021877e-04 -7,048175e-05
[1] "==== Covariâncias ===="
          [,1]      [,2]      [,3]
[1,] 1,0086483 0,9054011 0,8086425
[2,] 0,9054011 1,0029661 0,7056944
[3,] 0,8086425 0,7056944 1,0018611
[1] "==== Correlações ===="
          [,1]      [,2]      [,3]
[1,] 1,0000000 0,9001772 0,8044201
[2,] 0,9001772 1,0000000 0,7039954
[3,] 0,8044201 0,7039954 1,0000000
[1] "==== Covariâncias ===="
     [,1] [,2] [,3]
[1,]  1,0  0,9  0,8
[2,]  0,9  1,0  0,7
[3,]  0,8  0,7  1,0
In [5]:
#==============================================================================
# Utilizando a decomposicao de Cholesky
#==============================================================================
U     <- chol(SIGMA)        # decomposicao de SIGMA
y     <- z%*%t(U)           # gerando valores (note a transposicao)

print("==== Médias ====")
print(apply(x,2,"mean"))  # vetor (linhas) de medias amostral

print("==== Covariâncias ====")
print(var(x))             # matriz de covariancias amostral

print("==== Correlações ====")
print(cor(x))             # matriz de correlacoes amostral

print("==== Covariâncias ====")
print(SIGMA)
[1] "==== Médias ===="
[1]  3,225724e-06 -3,021877e-04 -7,048175e-05
[1] "==== Covariâncias ===="
          [,1]      [,2]      [,3]
[1,] 1,0086483 0,9054011 0,8086425
[2,] 0,9054011 1,0029661 0,7056944
[3,] 0,8086425 0,7056944 1,0018611
[1] "==== Correlações ===="
          [,1]      [,2]      [,3]
[1,] 1,0000000 0,9001772 0,8044201
[2,] 0,9001772 1,0000000 0,7039954
[3,] 0,8044201 0,7039954 1,0000000
[1] "==== Covariâncias ===="
     [,1] [,2] [,3]
[1,]  1,0  0,9  0,8
[2,]  0,9  1,0  0,7
[3,]  0,8  0,7  1,0