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.