[R-br] geração de variáveis aleatórias normais com estrutura de correlação

[Adriele Giaretta Biase]

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.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161017/26de2fe0/attachment.html>