[R-br] geração de variáveis aleatórias normais com estrutura de correlação
Éder Comunello
comunello.eder em gmail.com
Quarta Outubro 19 13:50:01 BRST 2016
Muito bom, Prof. Luiz Roberto!
Ficou muito didático dessa forma.
================================================
É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 18 de outubro de 2016 16:09, Luiz Roberto Martins Pinto <
luizroberto.uesc em gmail.com> escreveu:
> 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/20161019/38e483cb/attachment.html>
Mais detalhes sobre a lista de discussão R-br