| Ambos lados da revisão anteriorRevisão anteriorPróxima revisão | Revisão anterior |
| pessoais:eder [2011/09/26 08:01] – [section 5] eder | pessoais:eder [2018/07/04 19:10] (atual) – [section 5] eder |
|---|
| * {{:pessoais:reml_inla.r|Script}} Modelo seleção Genótipo ambiente via REML ML 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) | * {{:pessoais:linearregression.rnw|Script}} Regressão Linear - inferência via Mínimos quadrados, ML, REML, Gibbs, Metropolis, INLA, dclone ... (Em construção) |
| | * {{:pessoais:rjmcmc.r|RJMCMC}} Reversible Jump MCMC Regressão Linear |
| | ===== Artigos de Interesse ===== |
| | * {{http://www.sciencedirect.com/science/article/pii/S0022030278835905|Simulation of Examine Distributions of Estimators of Variances and Ratios of Variances}} |
| | * {{http://www.jstor.org/stable/3001853|Estimation of Variance and Covariance Components}} |
| | * {{http://www.sciencedirect.com/science/article/pii/S0022030291786013|C. R. Henderson: Contributions to Predicting Genetic Merit}} |
| | * {{http://www.sciencedirect.com/science/article/pii/S002203027584776X|Rapid Method for Computing the Inverse of a Relationship Matrix}} |
| | * {{http://www.jstor.org/stable/2527669?seq=18|The Estimation of Environmental and Genetic Trends from Records Subject to Culling}} |
| | * {{http://download.journals.elsevierhealth.com/pdfs/journals/0022-0302/PIIS002203027584776X.pdf|Rapid Method for Computing the Inverse of a Relationship Matrix}} |
| | * {{http://www.jstor.org/pss/2529339|A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values}} |
| | * {{http://www.jstor.org/stable/2529430?&Search=yes&searchText=%22C.+R.+Henderson%22&list=hide&searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&prevSearch=&item=3&ttl=373&returnArticleService=showFullText|Best Linear Unbiased Estimation and Prediction under a Selection Model}} |
| | * {{http://www.jstor.org/stable/2530609?&Search=yes&searchText=%22C.+R.+Henderson%22&list=hide&searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&prevSearch=&item=5&ttl=373&returnArticleService=showFullText|Variance-Covariance Matrix of Estimators of Variances in Unweighted Means ANOVA}} |
| | * {{http://www.jstor.org/stable/3001853?&Search=yes&searchText=%22C.+R.+Henderson%22&list=hide&searchUri=%2Faction%2FdoBasicSearch%3FQuery%3Dau%253A%2522C.%2BR.%2BHenderson%2522%26wc%3Don&prevSearch=&item=2&ttl=373&returnArticleService=showFullText| Estimation of Variance and Covariance Components}} |
| | |
| |
| ===== Disciplinas 2011/1 ===== | ===== Disciplinas 2011/1 ===== |
| * [[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: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 | * [[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 |
| | * [[http://www.leg.ufpr.br/doku.php/pessoais:eder:runicentro| Curso de Sofware R Analise Experimentos - UNICENTRO]] |
| ===== Códigos (Em construção) ===== | ===== Códigos (Em construção) ===== |
| <code R> | <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 | |
| |
| rj.modelo <- function(y,x,b0,b1,sigma,b01,b11,b21,sigma1,model,mu=mu,sd=sd,mu0=mu0, | |
| 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))} | |
| } | |
| |
| |
| | |
| 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])))) | |
| } | |
| |
| 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]]) | |
| ##------------------------------------------------------------------### | ##------------------------------------------------------------------### |
| ###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### |
