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:00Em 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/ USPContato: (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.