pessoais:eder - Wiki do LEG

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » pessoais » eder

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- pessoais:eder [2011/04/23 17:19]
eder [section 1]
+++ pessoais:eder [2011/09/26 08:01]
eder [section 5]
@@ Linha 1: / Linha 1: @@
 ===== Éder David Borges da Silva =====
-  * Graduando em Engenharia Agronômica - [[http://www.agronomia.ufpr.br|Eng Agronômica - UFPR]]
+  * Graduado em Engenharia Agronômica - [[http://www.agronomia.ufpr.br|Eng Agronômica - UFPR]]
   * **e-mail:** ederdbs@gmail.com / eder@leg.ufpr.br
@@ Linha 11: / Linha 11: @@
   *  Estatística Experimental
   *  Estatística Espacial
-  * [[http://www.leg.ufpr.br/doku.php/projetos:gem2|GEM²]] - em construção
+  * [[http://www.leg.ufpr.br/doku.php/projetos:gem2|GEM²]] Grupo de estudos em modelos mistos
+  * {{:pessoais:inlarrblup.r|GWS}} Seleção Genômica Ampla Via ML REML INLA
+  * {{:pessoais:reml_inla.r|Script}} Modelo seleção Genótipo ambiente via REML ML INLA
+  * {{:pessoais:linearregression.rnw|Script}} Regressão Linear - inferência via Mínimos quadrados, ML, REML, Gibbs, Metropolis, INLA, dclone ... (Em construção)
+===== Disciplinas 2011/1 =====
+  * [[http://www.leg.ufpr.br/doku.php/disciplinas:ce210-2010-02|CE-210: Inferência estatística II]]
+  * [[http://www.leg.ufpr.br/doku.php/disciplinas:ce718|CE-718: Métodos Computacionalmente Intensivos]]
 ===== Minicursos =====
    * [[http://www.leg.ufpr.br/ragronomia|Estatística Experimental com Software R]]
-   * [[http://www.leg.ufpr.br/doku.php/pessoais:eder:exptempo| Análise de Experimentos de longa duração]]
+   * [[http://www.leg.ufpr.br/doku.php/pessoais:eder:planejamentofito|Planejamento de experimento PG Produção Vegetal UFPR]]
+   * [[http://www.leg.ufpr.br/doku.php/pessoais:eder:exptempo| Análise de Experimentos de longa duração]] II Reunião Paranaense Ciência do Solo
+===== Códigos (Em construção) =====
+<code R>
+###-----------------------------------------------------------------###
+### Reversible jump MCMC
+### Modelo 1 y ~ N(b0+b1*x,sigma)
+### Modelo 1 y ~ N(b0+b1*x+b2*x²,sigma)
+#browseURL('http://www.icmc.usp.br/~ehlers/bayes/cap4.pdf')
+# pg 76
+require(MASS)#mvnorm()
+require(MCMCpack)#rinvgamma() dinvgamma()
+require(coda)#as.mcmc
+rm(list=ls())
+### Ajustar Prioris
+### conferir jacobiano
-===== Planejamento e Análise de experimentos =====
+rj.modelo <- function(y,x,b0,b1,sigma,b01,b11,b21,sigma1,model,mu=mu,sd=sd,mu0=mu0,
- == Material Geral do curso==
+                      mu02=mu02,V0=V0,V02=V02,v0=v0,tau0=tau0,v02=v02,tau02=tau02){
+    if (model == 1){
+      u <- rnorm(1, mu,sd)
+      b0_n <- b0
+      b1_n <- b1
+      sigma_n <- sigma
+      b01_n <- b0 * u
+      b11_n <-  b1 * u
+      b21_n <- u
+      sigma1_n <- sigma * (u^2)
+      }
+    if (model == 2){
+      u <- b21
+      b0_n <- b01 / u
+      b1_n <- b11 / u
+      sigma_n <-  sigma1 / (u^2)
+      b01_n <- b01
+      b11_n <-  b11
+      b21_n <- b21
+      sigma1_n <- sigma1
+      }
+    num <-  (sum(dnorm(y,b0_n+b1_n*x,sigma_n,log=TRUE))#+
+             #sum(dnorm(b0_n,mu0[1],V0[1,1],log=TRUE))+
+             #sum(dnorm(b1_n,mu0[1],V0[2,2],log=TRUE))+
+             #sum(log(dinvgamma(sigma_n,v0,tau0)))
+             ) * u^4
+    den <-  (sum(dnorm(y,b01_n+b11_n*x+b21_n*x^2,sigma1_n,log=TRUE))#+
+            # sum(dnorm(b01_n,mu02[1],V02[1,1],log=TRUE))+
+            # sum(dnorm(b11_n,mu02[2],V02[2,2],log=TRUE))+
+            # sum(dnorm(b21_n,mu02[3],V02[3,3],log=TRUE))+
+            # sum(log(dinvgamma(sigma1_n,v02,tau02)))
+             ) * dnorm(u,0,2)
+      u = runif(1, 0, 1)
+      if (model == 1) {
+        aceita = min(1, num/den)
+          if (u < aceita) {
+          model = 2
+          b0 <- b0_n
+          b1 <- b1_n
+          sigma <- sigma_n
+          }
+      }
+      if (model == 2){
+          aceita = min(1, den/num)
+            if (u < aceita) {
+              model = 1
+              b01 <- b01_n
+              b11 <- b11_n
+              b21 <- b21_n
+              sigma1 <- sigma1_n
+              }
+          }
+      if (model == 1){return(list(model = model,b0=b0,b1=b1,sigma=sigma))}
+      if (model == 2){return(list(model = model,b01=b01,b11=b11,b21=b21,sigma1=sigma1))}
+  }
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/CursoRFito.zip| Geral(Dados, Scripts...)]]
- == Revisão de Estatística Basica ==
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/RevisaoEstBasica.R| Script]]
+rjmcmc <- function(nI, x,y,burnIN,mu=mu,sd=sd) {
- == Planejamento de Experimentos ==
+  chain = matrix(NA, nrow = nI, ncol = 8)
-    *  [[http://www.leg.ufpr.br/~eder/CursoRFito/PlaneAnaExp.pdf| Apresentação]]
+  nv <- c(0,0)
- == Introdução ao R ==
+  chain[1,1:8] = c(1)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/IntroducaoR.R| Script]]
+  model = 1
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Introdu%e7%e3o%20ao%20Software%20R.pdf| Apresentação]]
+  n <- length(y)
- == Delineamento Inteiramente ao Acaso ==
+  ###----------------------------------------------------------###
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/DIC/exemplo%20dic.R|Script]]
+  ###MOdel 1
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/DIC/exemplo dic.xls|dados]]
+  X <- model.matrix(~x)
- == Delineamento Blocos ao Acaso ==
+  k<-ncol(X)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/DBA/exemplo%20dba.R|Script]]
+  #beta
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/DBA/exemplo dba.xls|dados]]
+  mu0<-rep(0,k)
- == Esquema Fatorial ==
+  V0<-100*diag(k)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Fatorial/Fatorial.R|Script]]
+  #sigma2
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Fatorial/fatorial2x3.txt|dados]]
+  v0<-3
- == Parcelas Sub-divididas ==
+  tau0<-100
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/SubDividida/Subdividida.R|Script]]
+  #Valores iniciais
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/SubDividida/subdividida.txt|dados]]
+  chain[1,3] <-sig2draw<- 3
- == Regressão ==
+  invV0 <- solve(V0)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao/Regre.r|Script]]
+  XtX <- crossprod(X,X)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao/Escoamento.txt|dados]]
+  Xty <- crossprod(X,y)
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao.pdf|Apresentação]]
+  invV0_mu0 <- invV0 %*% mu0
- == Script geral do curso ==
+  ###----------------------------------------------------------###
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Material%20SNW/Script.pdf| pdf]]
+  # Model 2
-    * [[http://www.leg.ufpr.br/~eder/CursoRFito/Script.R|Script]]
+  X2 <- cbind(1,x,x^2)
+  k2<-ncol(X2)
+  #beta
+  mu02<-rep(0,k2)
+  V02<-100*diag(k2)
+  #sigma2
+  v02<-3
+  tau02<-100
+  #Valores iniciais
+  chain[1,7] <- sig2draw2<- 3
+  invV02 <- solve(V02)
+  XtX2 <- crossprod(X2,X2)
+  Xty2 <- crossprod(X2,y)
+  invV0_mu02 <- invV02 %*% mu02
+  ###----------------------------------------------------------###
+  for (i in 2:nI) {
+    if (model == 1){
+         # Model 1
+         #beta
+         invsig2draw <- 1/sig2draw
+         V1<-solve(invV0+(invsig2draw) * XtX)
+         mu1<-V1 %*% (invV0_mu0 + (invsig2draw)* Xty)
+         chain[i,1:2]<-mvrnorm(n=1,mu1,V1)
+         # sigma
+         v1<-(n+2*v0)/2
+         yXb <- (y-X %*% chain[i,1:2])
+         tyXb <-t(yXb)
+         tau1<-(0.5)*(tyXb %*% yXb+2*tau0)
+         chain[i,3] <- sig2draw <- sqrt(rinvgamma(1,v1,tau1))
+      }
+    if (model == 2){
+         # Model 2
+         #beta
+         invsig2draw2 <- 1/sig2draw2
+         V12<-solve(invV02+(invsig2draw2) * XtX2)
+         mu12<-V12 %*% (invV0_mu02 + (invsig2draw2)* Xty2)
+         chain[i,4:6]<-mvrnorm(n=1,mu12,V12)
+         # sigma
+         v12<-(n+2*v02)/2
+         yXb2 <- (y-X2 %*% chain[i,4:6])
+         tyXb2 <-t(yXb2)
+         tau12<-(0.5)*(tyXb2 %*% yXb2+2*tau02)
+         chain[i,7] <- sig2draw2 <-  sqrt(rinvgamma(1,v12,tau12))
+        }
+    new <- rj.modelo(y,x,chain[i,1],chain[i,2],chain[i,3],chain[i,4],chain[i,5],chain[i,6],chain[i,7],model,mu=mu,sd=sd)
+    model  <-  new$model
+    if (model == 1) {
+            chain[i, 1] = new$b0
+            chain[i, 2] = new$b1
+            chain[i, 3] = new$sigma
+            nv[1] = nv[1] + 1
+            }
+    if (model == 2) {
+            chain[i, 4] = new$b01
+            chain[i, 5] = new$b11
+            chain[i, 6] = new$b21
+            chain[i, 7] = new$sigma1
+            nv[2] = nv[2] + 1
+    }
+}
+chain[,8] <- 1
+chain[is.na(chain[,1]),8] <- 2
+chain <- chain[- c(1:burnIN),]
+colnames(chain) <- c('b0_1','b1_1','sigma_1','b0_2','b1_2','b2_2','sigma_2','model')
+return(list(as.mcmc(na.omit(chain[,1:3])),
+            as.mcmc(na.omit(chain[,4:7])),
+            as.mcmc(na.omit(chain[,8]))))
+}
+x <- 1:10
+y <- 10+2*x^1+rnorm(x,0,5)
+plot(x,y)
+res <- rjmcmc(5000,x,y,1,mu=0,sd=100)
+lapply(res,summary)
+plot(res[[1]])
+summary(lm(y~1+I(x)))
+plot(res[[2]])
+summary(lm(y~1+I(x)+I(x^2)))
+plot(res[[3]])
+##------------------------------------------------------------------###
+###-----------------------------------------------------------------###
+### Regressão Beta
+### pacote oficial
+require(betareg)
+data("FoodExpenditure", package = "betareg")
+fe_beta <- betareg(I(food/income) ~ income + persons , data = FoodExpenditure)
+summary(fe_beta)
+###-----------------------------------------------------------------###
+### log vero da regressão beta com duas covariaveis,
+log.vero <- function(par,y,x1,x2){
+        mu <- exp((par[1] + par[2] * x1 + par[3] * x2))/(1+exp((par[1] + par[2] * x1 + par[3] * x2)))##logit^-1
+        ll  <- sum(dbeta(y, mu* par[4], (1-mu)*par[4],log = TRUE))
+        return(ll)
+}
+###-----------------------------------------------------------------###
+opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income,
+                                                        x1=FoodExpenditure$income,
+                                                        x2=FoodExpenditure$persons,
+                                                        hessian = TRUE, control=(list(fnscale=-1)))
+opt
+opt$par
+sqrt(-diag(solve(opt$hessian)))
+summary(fe_beta)
+###-----------------------------------------------------------------###
+log.veroP <- function(par,phi,y,x1,x2){
+        mu <- exp((par[1] + par[2] * x1 + par[3] * x2))/(1+exp((par[1] + par[2] * x1 + par[3] * x2)))##logit^-1
+        ll  <- sum(dbeta(y, mu* phi, (1-mu)*phi,log = TRUE))
+        return(ll)
+}
+opt <- grid.phi <- seq(20,60,l=150)
+con <- 1
+for (i in grid.phi){
+  opt[con] <- optim(c(B0=-0.5,B1=-0.51,B2=0.11),log.veroP,phi=i,y=FoodExpenditure$food/FoodExpenditure$income,
+                                                        x1=FoodExpenditure$income,
+                                                        x2=FoodExpenditure$persons,
+                                                        hessian = TRUE, control=(list(fnscale=-1)))$value
+  con <- con+1
+}
+plot(grid.phi,2*(max(opt)-opt),type='l')
+abline(h=3.84)
+</code>
+[[http://www.ime.usp.br/~sferrari/beta.pdf|Regressão beta]]

Diferenças

Navegação

Busca

Ferramentas

QR Code