Se você desejar entender o significado de correlação experimente:

n=10000000

cor_x1x2=-0.2

medx1=0.0067

varx1=0.0017

sdx1=varx1^0.5

medx2=0.1374

varx2=0.0024

sdx2=varx2^0.5

cov_x1x2=-cor_x1x2*(sdx1*sdx2)

sda=(varx1-cov_x1x2)^0.5

sdb=(varx2-cov_x1x2)^0.5

a=rnorm(n,medx1,sda);var(a)

b=rnorm(n,medx2,sdb);var(b)

c=rnorm(n,0,cov_x1x2^0.5);var(c)

x1=a+c

x2=b-c

var(x1)

var(x2)

cor(x1,x2)

Abraços,

Luiz Roberto.

Luiz Roberto Martins Pinto
Prof. Pleno/DCET/UESC

Laboratório de Estatística Computacional

Universidade Estadual de Santa Cruz

Ilhéus-Bahia-Brasil

luizroberto.uesc@gmail.com
skype: lrmpinto

http://lattes.cnpq.br/2732314327604831 

"The science exists because there are patterns. The patterns exist because God created them.

 The statistic exists to research the patterns that God created."

Em 18 de outubro de 2016 14:52, Éder Comunello via R-br <r-br@listas.c3sl.ufpr.br> escreveu:

Adriele, boa tarde!

Tentei adaptar seu código R com base no que entendi do script em SAS. Veja se pode ser útil...

### <code r>

library(MASS)
{
  set.seed(1357)
  n  <- 1000  # tamanho da amostra gerada (N)
  p  <- 2     # numero de variáveis a serem geradas (K)
  ME     <- rep(1, p)
  rho    <- -0.2 # correlacao negativa - 0.2
  sigma2 <- 1
  sigma  <- sigma2 * ((1-rho)*diag(p)+rho*matrix(1, p, p))
  SI  <- mvrnorm(n, ME, sigma)
}
apply(SI, 2, mean) # ~ 1
# [1] 1.037774 1.017846

cor(SI)            # ~ sigma
# [,1]       [,2]
# [1,]  1.0000000 -0.1973338
# [2,] -0.1973338  1.0000000

mu <- c(0.0067, 0.1374)
mean(0.0017*SI[,1]) # ~ 0.0017
mean(0.0024*SI[,2]) # ~ 0.0024

RES <- data.frame(X1=0.0067+0.0017*SI[,1], X2=0.1374+0.0024*SI[,2])

round(cbind(RES=colMeans(RES), INI=mu),4)
# RES    INI
# X1 0.0085 0.0067
# X2 0.1398 0.1374

cor(RES)
# X1         X2
# X1  1.0000000 -0.1973338
# X2 -0.1973338  1.0000000

### </code>

================================================
Éder Comunello
Researcher at Brazilian Agricultural Research Corporation (Embrapa)
DSc in Agricultural Systems Engineering (USP/Esalq)
MSc in Environ. Sciences (UEM), Agronomist (UEM)
---
Embrapa Agropecuária Oeste, Dourados, MS, Brazil |<O>|
================================================
GEO, -22.2752, -54.8182, 408m
UTC-04:00 / DST: UTC-03:00

Em 17 de outubro de 2016 16:05, Adriele Giaretta Biase via R-br <r-br@listas.c3sl.ufpr.br> escreveu:

Olá pessoal,

tenho uma dúvida com relação à geração de variáveis aleatórias normais multivariadas no R.

Eu gostaria de gerar duas variáveis aleatórias (x1, x2) usando distribuição normal multivariada com estrutura de correlação entre elas. A Estrutura da matriz de correlação foi estimada antes com base num banco de dados reais (com correlação fraca de - 0.2). Porém, essas variáveis não podem ser negativas, Ex. x1 tem média e variância, respectivamente, iguais a 0.0067 e 0.0017; x2 tem média e variância, respectivamente, iguais a 0.1374 e 0.0024. Eu gero as variáveis da seguinte forma:

Programa usado no R:

# Simulando os parâmetros com estrutura de correlação entre eles

n <- 1000  # tamanho da amostra gerada

p <- 2  # numero de variáveis a serem geradas

 

library(MASS)

# construíndo a matriz de correlação para usar nas simulações, baseadas na característica da amostra coletada

mu<-rep(0, times = p)

rho <-   - 0.2             # correlacao negativa - 0.2

sigma2 <- 1

Sigma <- sigma2 * ((1-rho)*diag(p)+rho*matrix(1, p, p))

X1 <- 0.0096  # Media da variavel x1

X2 <- 0.1203 # Media da variavel x2

media <- c(X1,X2)

y  <-  mvrnorm(n, media, Sigma)

cor(y)

apply(y, 2, mean)

y

   [1,]  0.309910452  1.0642521083

   [2,] -0.251583312  1.8909032058

   [3,]  1.330362012 -1.1239501814

   [4,] -0.793399464 -1.5433056284

   [5,]  2.165144843 -0.2645184534

   [6,]  0.532777085  1.0910864562

   [7,] -1.612135390  2.0489354648

   [8,] -0.430529913  1.0312062602

...

 

Quando gero as simulações para x1 e x2 usando a função mvrnorm, o resultado me retorna alguns valores negativos para as variáveis, isso não poderia ocorrer. Teria alguma outra função em que eu possa fornecer a variância de cada uma dessas variáveis, além da estrutura de correlação? Como poderia contornar essa situação, usando o R?

 

Obs: No SAS, a seguinte programação funciona perfeitamente:

proc iml;

wrksize=100000;

K=2;

N=1000; /* tamanho da amostra */

M={0 0};

S={ 1       -0.5283,

   -0.5283       1 };

X=shape(0,K,N);

ME=0; SI=1;

DO I=1 TO K;

  DO J=1 TO N;

    if I>1

      then

        do;

         ME=M[I]+(S[1:I-1,I])`*(inv(S[1:I-1,1:I-1])*(X[1:I-1,J]-M[1:I-1]));

         SI=S[I,I]-(S[1:I-1,I])`*(inv(S[1:I-1,1:I-1])*(S[1:I-1,I]));

        end;

    X[I,J]=ME+NORMAL(0)*SQRT(SI);

 

  END;

 

END;

Z=t(X);

 

varnames='X1':'X2';

create NOVO from Z[colname=varnames];

append from Z;

quit;

 

data MEXT; set NOVO;

options ps=66 ls=75;

   X1=0.0067+0.0017*X1; /* normal */

   X2=0.1374+0.0024*X2; /* normal */

 

proc corr data=MEXT;

var X1 X2;

run;

--

Adriele Giaretta Biase.

Mestre em  Estatística e Experimentação Agropecuária - UFLA. 
Doutora em Estatística e Experimentação Agronômica - ESALQ/ USP

Contato: (19) 98861-0619.

_______________________________________________
R-br mailing list
R-br@listas.c3sl.ufpr.br
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.

_______________________________________________
R-br mailing list
R-br@listas.c3sl.ufpr.br
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.