Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
pessoais:eder [2011/06/01 22:31] eder [Codigos] |
pessoais:eder [2011/09/25 20:50] eder [Códigos] |
||
---|---|---|---|
Linha 12: | Linha 12: | ||
* Estatística Espacial | * Estatística Espacial | ||
* [[http://www.leg.ufpr.br/doku.php/projetos:gem2|GEM²]] Grupo de estudos em modelos mistos | * [[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 ===== | ===== Minicursos ===== | ||
* [[http://www.leg.ufpr.br/ragronomia|Estatística Experimental com Software R]] | * [[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]] |
- | ===== Codigos ===== | + | * [[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> | <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 | ||
+ | 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))} | ||
+ | } | ||
+ | |||
| | ||
- | buf(100000) | + | 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 | ||
+ | } | ||
+ | } | ||
+ | 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])))) | ||
+ | } | ||
- | ### MOnte carlo | + | x <- 1:10 |
- | ### inversão de p | + | y <- 10+2*x^1+rnorm(x,0,5) |
- | ################################################################################ | + | plot(x,y) |
- | ### Regressão beta | + | res <- rjmcmc(5000,x,y,1,mu=0,sd=100) |
- | rm(list=ls()) | + | lapply(res,summary) |
- | require(betareg) | + | 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 | ### pacote oficial | ||
+ | require(betareg) | ||
data("FoodExpenditure", package = "betareg") | data("FoodExpenditure", package = "betareg") | ||
fe_beta <- betareg(I(food/income) ~ income + persons , data = FoodExpenditure) | fe_beta <- betareg(I(food/income) ~ income + persons , data = FoodExpenditure) | ||
summary(fe_beta) | summary(fe_beta) | ||
- | ###----------------------------------------------------------### | + | ###-----------------------------------------------------------------### |
- | ### log vero | + | ### log vero da regressão beta com duas covariaveis, |
- | log.vero <- function(B0,B1,B2,phi,y,x1,x2){ | + | log.vero <- function(par,y,x1,x2){ |
- | mu <- exp((B0 + B1 * x1 + B2 * x2))/(1+exp((B0 + B1 * x1 + B2 * x2)))##logit^-1 | + | 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)) | + | ll <- sum(dbeta(y, mu* par[4], (1-mu)*par[4],log = TRUE)) |
return(ll) | return(ll) | ||
} | } | ||
- | ###----------------------------------------------------------### | + | |
- | log.vero(-0.62,-0.12,0.11,35,y=FoodExpenditure$food/FoodExpenditure$income, | + | ###-----------------------------------------------------------------### |
- | x1=FoodExpenditure$income, | + | opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income, |
- | x2=FoodExpenditure$persons) | + | |
- | ###----------------------------------------------------------### | + | |
- | ### B0 B1 | + | |
- | par.vals <- expand.grid(B0=seq(0,2,l=100),B1=seq(-1,1,l=100)) | + | |
- | logL <- apply(as.matrix(par.vals),1,log.vero,B2=0.11,phi=35,y=FoodExpenditure$food/FoodExpenditure$income, | + | |
- | x1=FoodExpenditure$income, | + | |
- | x2=FoodExpenditure$persons) | + | |
- | contour(unique(par.vals$B0),unique(par.vals$B1),matrix(logL,ncol=100)) | + | |
- | ###----------------------------------------------------------### | + | |
- | opt <- optim(c(-0.5,-0.12,0.11,35),logvero,y=FoodExpenditure$food/FoodExpenditure$income, | + | |
x1=FoodExpenditure$income, | x1=FoodExpenditure$income, | ||
x2=FoodExpenditure$persons, | x2=FoodExpenditure$persons, | ||
hessian = TRUE, control=(list(fnscale=-1))) | 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> | </code> | ||
+ | |||
[[http://www.ime.usp.br/~sferrari/beta.pdf|Regressão beta]] | [[http://www.ime.usp.br/~sferrari/beta.pdf|Regressão beta]] | ||
- | ===== Planejamento e Análise de experimentos ===== | ||
- | == Material Geral do curso== | ||
- | |||
- | * [[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]] | ||
- | == Planejamento de Experimentos == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/PlaneAnaExp.pdf| Apresentação]] | ||
- | == Introdução ao R == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/IntroducaoR.R| Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Introdu%e7%e3o%20ao%20Software%20R.pdf| Apresentação]] | ||
- | == Delineamento Inteiramente ao Acaso == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/DIC/exemplo%20dic.R|Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/DIC/exemplo dic.xls|dados]] | ||
- | == Delineamento Blocos ao Acaso == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/DBA/exemplo%20dba.R|Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/DBA/exemplo dba.xls|dados]] | ||
- | == Esquema Fatorial == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Fatorial/Fatorial.R|Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Fatorial/fatorial2x3.txt|dados]] | ||
- | == Parcelas Sub-divididas == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/SubDividida/Subdividida.R|Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/SubDividida/subdividida.txt|dados]] | ||
- | == Regressão == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao/Regre.r|Script]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao/Escoamento.txt|dados]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Regressao.pdf|Apresentação]] | ||
- | == Script geral do curso == | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Material%20SNW/Script.pdf| pdf]] | ||
- | * [[http://www.leg.ufpr.br/~eder/CursoRFito/Script.R|Script]] | ||
- | |||