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.
