[R-br] geração de variáveis aleatórias normais com estrutura de correlação
Luiz Roberto Martins Pinto
luizroberto.uesc em gmail.com
Terça Outubro 18 17:09:12 BRST 2016
Se você desejar entender o significado de correlação experimente:
n=10000000
cor_x1x2=-0.2
medx1=0.0067
varx1=0.0017
sdx1=varx1^0.5
medx2=0.1374
varx2=0.0024
sdx2=varx2^0.5
cov_x1x2=-cor_x1x2*(sdx1*sdx2)
sda=(varx1-cov_x1x2)^0.5
sdb=(varx2-cov_x1x2)^0.5
a=rnorm(n,medx1,sda);var(a)
b=rnorm(n,medx2,sdb);var(b)
c=rnorm(n,0,cov_x1x2^0.5);var(c)
x1=a+c
x2=b-c
var(x1)
var(x2)
cor(x1,x2)
Abraços,
Luiz Roberto.
Luiz Roberto Martins Pinto
Prof. Pleno/DCET/UESC
Laboratório de Estatística Computacional
Universidade Estadual de Santa Cruz
Ilhéus-Bahia-Brasil
luizroberto.uesc em gmail.com
skype: lrmpinto
http://lattes.cnpq.br/2732314327604831
"*The s**cience exists because there are patterns. The patterns exist
because God created them*.
* The statistic exists to research the patterns that God created.*"
Em 18 de outubro de 2016 14:52, Éder Comunello via R-br <
r-br em listas.c3sl.ufpr.br> escreveu:
> 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.
>>
>
>
> _______________________________________________
> 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/d54ee426/attachment-0001.html>
Mais detalhes sobre a lista de discussão R-br