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

[Éder Comunello]

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:00




Em 17 de outubro de 2016 16:05, Adriele Giaretta Biase via R-br <
r-br em 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/ 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/8311c3ef/attachment.html>