mec {geoComp}R Documentation

Usage

mec(dados.comp, metodo = "L-BFGS-B", otimizador = "optim", print.pars = FALSE, alpha = 0.975)

Arguments

dados.comp
metodo
otimizador
print.pars
alpha

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--    or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function(dados.comp, metodo="L-BFGS-B",otimizador="optim",print.pars=FALSE, alpha=0.975){
if(class(dados.comp) != "geoComp") stop("A classe do objeto deve ser geoComp")
y1 <- dados.comp[[1]]$Y1
y2 <- dados.comp[[1]]$Y2
coords <- dados.comp[[3]]
X <- cbind(rep(1:0,length=length(y1)*2),rep(0:1,length=length(y1)*2))
  
## 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)

### Pegando o Y
y <- dados.comp[[2]][[1]]

## Reparametrizando
eta <- s2/s1
nu1 <- tau1/s1
nu2 <- tau2/s1
theta1 <- c(eta,nu1,nu2,phi,rho)

## Calculando a logvero para o chute inicial
ll <- log.vero(theta1,dados.comp)

## Otimizando
if(metodo== "L-BFGS-B"){
   estim <- optim(theta1,log.vero,dados.comp = dados.comp, hessian=TRUE,method=metodo,print.pars=print.pars,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,dados.comp = dados.comp,hessian=TRUE,method=metodo,print.pars=print.pars)
   ifelse(estim$convergence[1] == 0, print("O algoritmo convergiu"),print("Existe problemas na convergencia"))}

## Voltando aos valores iniciais
thetaest <- estim$par
Vest <- monta.V(thetaest,dados.comp)
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)))
## Intervalo de confianca pelo metodo Delta
ic.min <- estima - errospd
ic.max <- estima + errospd
retorna <- list()
retorna[1][[1]] <- data.frame(para,estima,errospd,ic.min,ic.max)
names(retorna[1][[1]]) <- c("Parametros","Estimativas","Erro Padrao","LI.Delta","LS.Delta")
retorna[2][[1]] <-c(ll,llest) 
names(retorna[2][[1]]) <- c("LogLik Start", "LogLik Optim")
return(retorna)
  }

[Package geoComp version 0.1-0 Index]