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/06/05 19:50]
eder
+++ pessoais:eder [2011/09/25 20:50]
eder [Códigos]
@@ Linha 12: / Linha 12: @@
   *  Estatística Espacial
   * [[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
+===== 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]]
-===== Códigos =====
+   * [[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>
-###buf
+###-----------------------------------------------------------------###
-buf <- function(n){
+### Reversible jump MCMC
-  ttt <- NULL
+### Modelo 1 y ~ N(b0+b1*x,sigma)
-  ttt[1] <- 0
+### Modelo 1 y ~ N(b0+b1*x+b2*x²,sigma)
-  x <- runif(n)
+#browseURL('http://www.icmc.usp.br/~ehlers/bayes/cap4.pdf')
-  th <- runif(n,0,pi)
+# pg 76
-  st <- sin(th)
+require(MASS)#mvnorm()
-  for ( i in 1:n){
+require(MCMCpack)#rinvgamma() dinvgamma()
-    if(st[i]>x[i]){
+require(coda)#as.mcmc
-      ttt[i+1]  <- ttt[i]+1
+rm(list=ls())
-    }
+### Ajustar Prioris
-    else {
+### conferir jacobiano
-      ttt[i+1] <- ttt[i]
-    }}
+rj.modelo <- function(y,x,b0,b1,sigma,b01,b11,b21,sigma1,model,mu=mu,sd=sd,mu0=mu0,
-    if (ttt[n+1]>0){
+                      mu02=mu02,V0=V0,V02=V02,v0=v0,tau0=tau0,v02=v02,tau02=tau02){
-      plot((0:n)[ttt>0],2*(0:n)[ttt>0]/ttt[ttt>0],type='l',xlab='numero simulação',ylab='pi')
+    if (model == 1){
-    }
+      u <- rnorm(1, mu,sd)
-    else{print('no sucesso')}
+      b0_n <- b0
-    abline(pi,0)
+      b1_n <- b1
-    }
+      sigma_n <- sigma
+      b01_n <- b0 * u
-  buf(100000)
+      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(y,mu0[1],V0[1,1],log=TRUE))+
+             #sum(dnorm(y,mu0[1],V0[2,2],log=TRUE))+
+            # sum(log(dinvgamma(y,v0,tau0)))
+             ) * u^4
+    den <-  (sum(dnorm(y,b01_n+b11_n*x+b21_n*x^2,sigma1_n,log=TRUE))#+
+            # sum(dnorm(y,mu02[1],V02[1,1],log=TRUE))+
+            # sum(dnorm(y,mu02[2],V02[2,2],log=TRUE))+
+            # sum(dnorm(y,mu02[3],V02[3,3],log=TRUE))+
+            # sum(log(dinvgamma(y,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))}
+  }
-### MOnte carlo
-## Calcula a área via simulação de monte carlo
-## args: r= raio, s vetor com numero de simulação, plotS plotar a simulação
-MCcirculo<-function(r,s,plotS=TRUE){
-ns<-area<-s
-r<-r
-con <- 1
-for (j in ns) {
-#pontos aleatorios
-	x<-runif(j, min=-r, max=r)
-	y<-runif(j, min=-r, max=r)
-	ponto<-cbind(x,y)
-  cont <- sum(apply(ponto,1,function(x){sqrt(sum(x^2))})<r)
-#plotando Simulação
-  if(plotS==TRUE){
-	plot(x,y,col="red",type="p",asp=1,lwd=1,xlim=c(-r,r),ylim=c(-r,r), main="Simulação Monte Carlo",sub=j)
-	ang <- seq(0, 2*pi, length = 100)
-	xx <- r * cos(ang);yy <- r * sin(ang)
-	polygon(xx, yy,border = "dark blue",lwd=2)
-  }
-#Calculo de Area
-	area[con]<-(cont/j)*(r^2)*4
-  cat(paste(round(area[con],6),j,'\n'))
-  con <- con+1
-}
-	plot(ns,area,main="Simulação Monte Carlo",xlab='Número da amostra',ylab='Area')
-  abline(h=pi*r^2,col='red',lwd=2)
+rjmcmc <- function(nI, x,y,burnIN,mu=mu,sd=sd) {
+  chain = matrix(NA, nrow = nI, ncol = 8)
+  nv <- c(0,0)
+  chain[1,1:8] = c(1)
+  model = 1
+  n <- length(y)
+  ###----------------------------------------------------------###
+  ###MOdel 1
+  X <- model.matrix(~x)
+  k<-ncol(X)
+  #beta
+  mu0<-rep(0,k)
+  V0<-100*diag(k)
+  #sigma2
+  v0<-3
+  tau0<-100
+  #Valores iniciais
+  chain[1,3] <-sig2draw<- 3
+  invV0 <- solve(V0)
+  XtX <- crossprod(X,X)
+  Xty <- crossprod(X,y)
+  invV0_mu0 <- invV0 %*% mu0
+  ###----------------------------------------------------------###
+  # Model 2
+  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
+    }
 }
-MCcirculo(1,seq(5,5000,by=1000),plotS=FALSE)
+chain[,8] <- 1
-### inversão de p
+chain[is.na(chain[,1]),8] <- 2
-### Inversão de Probabilidade
+chain <- chain[- c(1:burnIN),]
-NS <- 10000
+colnames(chain) <- c('b0_1','b1_1','sigma_1','b0_2','b1_2','b2_2','sigma_2','model')
-U <- runif(NS)
+return(list(as.mcmc(na.omit(chain[,1:3])),
-X <- - log(U)
+            as.mcmc(na.omit(chain[,4:7])),
-Y <- rexp(NS)
+            as.mcmc(na.omit(chain[,8]))))
-par(mfrow=c(1,3))
+}
-hist(U,freq=FALSE,main='Uniforme',col='lightblue')
-lines(density(U),col='red',lwd=2)
+x <- 1:10
-hist(X,freq=FALSE,main='Expoencial via uniforme',col='lightblue')
+y <- 10+2*x^1+rnorm(x,0,5)
-lines(density(X),col='red',lwd=2)
+plot(x,y)
-lines(curve(dexp(x,1),min(X),max(X),add=TRUE),col='blue',lwd=2)
+res <- rjmcmc(5000,x,y,1,mu=0,sd=100)
-hist(Y,freq=FALSE,main='Expoencial do R',col='lightblue')
+lapply(res,summary)
-lines(density(Y),col='red',lwd=2)
+plot(res[[1]])
-lines(curve(dexp(x,1),min(Y),max(Y),add=TRUE),col='blue',lwd=2)
+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
@@ Linha 94: / Linha 208: @@
 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){
@@ Linha 101: / Linha 215: @@
         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,
@@ Linha 111: / Linha 225: @@
 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