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 | ||
pessoais:eder [2011/09/26 08:01] eder [section 5] |
pessoais:eder [2018/07/04 19:10] (atual) eder [section 5] |
||
---|---|---|---|
Linha 15: | Linha 15: | ||
* {{: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 ===== | ||
Linha 25: | Linha 38: | ||
* [[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]]) | ||
##------------------------------------------------------------------### | ##------------------------------------------------------------------### | ||
###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### |