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

[Elias T. Krainski]

O erro está na definicao de Sigma, você fez os cálculos com variâncias 
marginais igual a 1 e segundo o "enunciado" essas variâncias são bem 
menores.


On 17/10/2016 17:05, Adriele Giaretta Biase via R-br wrote:
>
> 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*;
>
> DOI=*1*TOK;
>
> DOJ=*1*TON;
>
> ifI>*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';
>
> createNOVO fromZ[colname=varnames];
>
> appendfromZ;
>
> *quit*;
>
> *data*MEXT; setNOVO;
>
> optionsps=*66*ls=*75*;
>
>    X1=*0.0067*+*0.0017**X1; /* normal */
>
> X2=*0.1374*+*0.0024**X2; /* normal */
>
> *proc**corr*data=MEXT;
>
> varX1 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 em 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.

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161018/3701db15/attachment.html>