Diferenças

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

--- pessoais:eder [2011/09/25 20:57] – [Área de Interesse] eder
+++ pessoais:eder [2018/07/04 19:10] (atual) – [section 5] eder
@@ Linha 15: / Linha 15: @@
   * {{: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: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 =====
@@ 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: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) =====
 <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(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))}
-  }
-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]])
 ##------------------------------------------------------------------###
 ###-----------------------------------------------------------------###

Diferenças

Navegação

Busca

Ferramentas

QR Code