infbayes <- function(dados,tune,queima,salto,nsim,verbose=TRUE,mediahiper){ y1 <- dados[[1]]$Y1 y2 <- dados[[1]]$Y2 coords <- dados[[3]] mu1 <- mean(y1) mu2 <- mean(y2) var_y1 <- var(y1) s1 <- var_y1/2 tau1 <- s1 var_y2 <- var(y2) s2 <- var_y2/2 tau2 <- s2 dim <- range(dist(coords)) phi <- dim[1] + 0.2*(dim[2]-dim[1]) rho <- cor(y1,y2) eta <- s2/s1 nu1 <- tau1/s1 nu2 <- tau2/s1 theta1 <- c(eta,nu1,nu2,phi,rho) if(verbose==TRUE)print("Otimizando a log-verossimilhanca") inicial.optim <- optim(theta1, logarit.vero, dados.comp=dados, method="L-BFGS-B", lower=c(1e-32,1e-32,1e-32,1e-32,-1), upper=c(Inf,Inf,Inf,Inf,1),hessian=TRUE) print(inicial.optim$par) inv.hessiana <-solve(inicial.optim$hessian) sd.5 <- sqrt(c(inv.hessiana[1,1],inv.hessiana[2,2],inv.hessiana[3,3], inv.hessiana[4,4],inv.hessiana[5,5])) tune <- tune diag.tune <- diag(c(tune)) mc.V <- diag.tune %*% inv.hessiana %*% diag.tune if(verbose==TRUE)print("Rodando o MCMC") theta1.inicial <- inicial.optim$par if(length(mediahiper)==0){ mh <- MCMCmetrop1R(posteriori,theta1.inicial, mediahiper=c(theta1.inicial[1:4],sd.5[1:4]), dados.comp=dados,burnin=queima,mcmc=nsim,thin=salto, V=mc.V) } if(length(mediahiper)==8){ mh1 <- MCMCmetrop1R(posteriori,theta1.inicial, mediahiper=c(mediahiper[1:4],mediahiper[5:8]), dados.comp=dados,burnin=queima,mcmc=nsim,thin=salto, V=mc.V) } #plot(mh) saida.mh <<- mh n <- length(dados[[2]]$Y) p <- length(dados[[1]]) D <- cbind(rep(1:0,length=n),rep(0:1,length=n)) sigma <- c() tamanho <- dim(mh)[1] sigma.media <- matrix(ncol=3,nrow=tamanho) for(i in 1:tamanho){ V <- monta.V(mh[i,],dados) # beta nao eficiente # inv.V <- solve(V) # beta <- solve(t(D)%*%inv.V%*%D)%*%(t(D)%*%inv.V%*%dados[[2]]$Y) # beta eficiente Rbayes <- chol(V) back.bayes1 <- backsolve(Rbayes,D,upper.tri=TRUE,transpose=TRUE) back.bayes2 <- backsolve(Rbayes,dados[[2]]$Y,upper.tri=TRUE,transpose=TRUE) inv.tDinvVD <- solve(crossprod(back.bayes1,back.bayes1)) beta <- inv.tDinvVD%*%crossprod(back.bayes1,back.bayes2) # s2 nao eficiente: # s2 <- 1/(n-p)*t(dados[[2]]$Y-D%*%beta)%*%inv.V%*%(dados[[2]]$Y-D%*%beta) # s2 eficiente: desvio.bayes <- dados[[2]]$Y-D%*%beta back.desvio.bayes <- backsolve(Rbayes,desvio.bayes, upper.tri=TRUE,transpose=TRUE) s2 <- 1/(n-p)*crossprod(back.desvio.bayes,back.desvio.bayes) sigma[i] <- rinvchisq(1,n-p,s2) # var.beta nao eficiente: # var.beta <- as.numeric(sigma[i])*solve(t(D)%*%inv.V%*%D) # var.beta eficiente: var.beta <- as.numeric(sigma[i])*inv.tDinvVD sigma.media[i,] <- c(mvrnorm(1,mu=beta,Sigma=var.beta),sqrt(sigma[i])) print(i) } retorna <- list() retorna[[1]] <- data.frame(sigma.media) names(retorna[[1]]) <- c('beta1','beta2','sigma1') retorna[[2]] <- data.frame(mh) names(retorna[[2]]) <- c('eta','nu1','nu2','phi','rho') retorna[[3]] <- data.frame(inicial.optim$par) return(retorna) }