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 Ambos lados da revisão seguinte | ||
|
pessoais:eder [2011/09/26 08:01] eder [section 5] |
pessoais:eder [2011/09/26 20:55] eder [section 5] |
||
|---|---|---|---|
| Linha 27: | Linha 27: | ||
| ===== 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]]) | ||
| ##------------------------------------------------------------------### | ##------------------------------------------------------------------### | ||
| ###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### | ||
