# Ref:mgec.R ### Funcao para calcular o modelo espacial composicional########## mec <- function(y1,y2,X,coords,metodo="BFGS",otimizador="optim",alpha=0.975){ ## Calculando os valores iniciais ## 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)) # Um "chute" para phi phi <- dim[1] + 0.2*(dim[2]-dim[1]) rho <- cor(y1,y2) ### Organizando os dados seq11 <- seq(1,length(y1)*2,by=2) seq22 <- seq(2,length(y2)*2,by=2) y <- c() y[seq11] <- y1 y[seq22] <- y2 dados[1][[1]] <- y dados[2][[1]] <- coords ## Reparametrizando eta <- s2/s1 nu1 <- tau1/s1 nu2 <- tau2/s1 theta1 <- c(eta,nu1,nu2,phi,rho) ## Chamando ll da funcao log.vero ll <- log.vero(theta1) ## Otimizando if(otimizador == "optim"){ if(metodo== "L-BFGS-B"){ estim <- optim(theta1,log.vero,hessian=TRUE,method=metodo,print.pars=T,lower=c(1e-32,1e-32,1e-32,1e-32,-1),upper=c(Inf,Inf,Inf,Inf,1)) ifelse(estim$convergence[1] == 0, print("O algoritmo convergiu"),print("Existe problemas na convergencia")) } if(metodo=="Nelder-Mead" | metodo=="CG" | metodo=="BFGS"){ estim <- optim(theta1,log.vero,hessian=TRUE,method=metodo,print.pars=T) ifelse(estim$convergence[1] == 0, print("O algoritmo convergiu"),print("Existe problemas na convergencia")) } } if(otimizador == "nlminb"){ #fdhess<- fdHess(c(12.3, 2.34), function(x) x[1]*(1-exp(-0.4*x[2]))) #fdhess<- fdHess(theta1, function(theta1) log.vero(theta1)) estim <- nlminb(theta1,log.vero,hessian=fdhess,lower=c(1e-32,1e-32,1e-32,1e-32,-1),upper=c(Inf,Inf,Inf,Inf,1),print.pars=T) ifelse(estim$convergence[1] == 0, print("O algoritmo convergiu"),print("Existe problemas na convergencia")) } ## Voltando aos valores iniciais X <- cbind(rep(1:0,length=length(y1)*2),rep(0:1,length=length(y1)*2)) thetaest <- estim$par Vest <- monta.V(thetaest,dados) ldetVest <- determinant(Vest,log=TRUE)$modulus[1] muest<- drop(solve(crossprod(X,solve(Vest,X)))%*%crossprod(X,solve(Vest,y))) Qeest <- drop(crossprod(y,solve(Vest,y))-2*crossprod(y,solve(Vest,X%*%muest))+crossprod(muest,crossprod(X,solve(Vest,X%*%muest)))) n <- length(y) s1est <- sqrt(Qeest/n) llest <- drop(-0.5*(n*log(2*pi)+n*log(s1est^2)+ ldetVest + (1/(s1est^2))*Qeest)) s2est <- thetaest[1]*s1est tau1est <- thetaest[2]*s1est tau2est <- thetaest[3]*s1est ## Calculando a variāncia para mu e s1: 1/Io(mu) - 1/Io(s1) V.g.mu <- (s1est)^2*solve(crossprod(X,solve(Vest,X))) V.g.s1 <- (s1est)^3/(3*Qeest-n*s1est) ## Calculando o gradiente g<- matrix(c(s1est,0,0,0,0,0,s1est,0,0,0,0,0,s1est,0,0,0,0,0,1,0,0,0,0,0,1),nr=5) ## Calculando a variancia para s2,tau1,tau2,phi e rho - usa hessiana V.g <- crossprod(g,diag(1/diag(estim$hessian))%*%g) ## Preparando a saida para <- c("mu1","mu2","s1","s2","tau1","tau2","phi","rho") estima <- c(muest[1],muest[2],s1est,s2est,tau1est,tau2est,thetaest[4],thetaest[5]) errospd <- qnorm(alpha)*sqrt(c(V.g.mu[1],V.g.mu[2],V.g.s1,diag(V.g))) ic.min <- estima - errospd ic.max <- estima + errospd retorna <- list() retorna[1][[1]] <- data.frame(para,estima,errospd,ic.min,ic.max) retorna[2][[1]] <-c(ll,llest) return(retorna) }