# assis load("P20eea2.Rdata") P20eea2<-P20eea2[,c(-7,-8)] head(P20eea2) # estacional load("P20eec2.Rdata") P20eec2<-P20eec2[,c(-7,-8)] head(P20eec2) # cb load("P20pecb2.Rdata") head(P20pecb2) # ic load("P20peic2.Rdata") P20peic2<-P20peic2[,c(-7)] head(P20peic2) # 2 a 2 Paic<-rbind(P20eea2,P20peic2) Pccb<-rbind(P20eec2,P20pecb2) # todos Pt<-rbind(P20eea2,P20peic2,P20eec2,P20pecb2) require(geoR) require(RandomFields) a<-as.geodata(P20eea2,coords.col=2:3,data.col=4) c<-as.geodata(P20eec2,coords.col=2:3,data.col=4) cb<-as.geodata(P20pecb2,coords.col=2:3,data.col=4) ic<-as.geodata(P20peic2,coords.col=2:3,data.col=4) aic<-as.geodata(Paic,coords.col=2:3,data.col=4) ccb<-as.geodata(Pccb,coords.col=2:3,data.col=4) t<-as.geodata(Pt,coords.col=2:3,data.col=4) likfit(a, ini=c(sa, phia),nug=ta) likfit(c, ini=c(sc, phic),nug=tc) likfit(cb, ini=c(scb, phicb),nug=tcb) likfit(ic, ini=c(sic, phiic),nug=tic) #parT<-(mu, s2, tau2, phi ) parTa <- c(22, 5, 13, 76) parTc <- c(25, 115, 0, 4) parTcb <- c(27, 170, 157, 33) parTic <- c(29, 8, 68, 109) # Medias, log(sqrt(s2)), tau2, phi #parT<-(mu, s2, phi, tau2 ) # 1 a 1 # assis ma<-mean(a$data) s2a<- parTa[2] sa<-log(sqrt(s2a)) ta<-parTa[3] phia<-parTa[4] ma; s2a; sa; ta; phia # estacional mc<-mean(c$data) s2c<-parTc[2] sc<-log(sqrt(s2c)) tc<-parTc[3] phic<-parTc[4] mc; s2c; sc; tc; phic # Ombrofila mcb<-mean(cb$data) s2cb<- parTcb[2] scb<-log(sqrt(s2cb)) tcb<-parTcb[3] phicb<-parTcb[4] mcb; s2cb; tcb; phicb # Restinga mic<-mean(ic$data) s2ic<-parTic[2] sic<-log(sqrt(s2ic)) tic<-parTic[3] phiic<-parTic[4] mic; s2ic; tic; phiic # 2 a 2 # assis + restinga maic<- mean(aic$data) saic<- log(sqrt(s2a+s2ic/2)) taic<-(ta+tic)/2 phiaic<- (phia+phiic)/2 maic; saic; taic; phiaic # Ombrofila + estacional mccb<-mean(ccb$data) sccb<- log(sqrt(s2c+s2cb/2)) tccb<-(tc+tcb)/2 phiccb<- (phic+phicb)/2 mccb; sccb; tccb; phiccb # todos # assis + estacional + ombrofila + restinga mt<- mean(t$data) st<- log(sqrt(s2a+s2c+s2cb+s2ic/4)) tt<-(ta+tc+tcb+tic)/4 phit<- (phia+phic+phicb+phiic)/4 mt;st;tt;phit ### Parametros - importante # tau é igual tausq # log(sqrt(tau2))=par[1] # sqrt(tau2)=exp(par[1]) # tau2 = exp(p[1])^2 nloglik <- function(par, dados, dmin, print.pars=TRUE) { tau <- exp(par[1])^2 s2 <- exp(par[2])^2 phi <- exp(par[3]) u=as.matrix(dist(a$coords)) n <- length(dados$data) if (phi<0.001*dmin | s2<0.001*var(dados$data)) { ldet <- n*log(tau + s2) edet <-exp(ldet) s2 <- 0 phi <- 0 if(ldet=="Inf") return(ldet=.Machine$double.xmax^.5) if(edet=="Inf") return(edet=.Machine$double.xmax^.5) # ldet <- .Machine$double.xmax^.5 ### logL <- sum(dnorm(dados$data, mean(dados$data), tau+s2, log=TRUE)) Sig <- diag(tau+s2, n) } else { Sig <- diag(tau, n) + s2*exp(-u/phi) ldet <- log(det(Sig)) edet <-exp(ldet) if(ldet=="Inf") return(ldet=.Machine$double.xmax^.5) if(edet=="Inf") return(edet=.Machine$double.xmax^.5) # ldet <- .Machine$double.xmax^.5 } w <- solve(Sig) x <- matrix(1, nrow=nrow(dados$coords)) mu <- solve(t(x)%*%w%*%x, t(x)%*%w%*%dados$data) logL <- (-n/2)*log(2*pi)- (1/2)*ldet - ((1/2)*(sum((dados$data-mu) * solve(Sig, dados$data-mu)))) if(print.pars) print(c(exp(ldet),ldet,logL, par)) return(-logL) } mufun <- function(par, dados) { n <- nrow(dados$coords) x <- matrix(1, nrow=n) u <- as.matrix(dist(dados$coords)) Sig <- diag(par[1], n) + par[2]*exp(-u/par[3]) w <- solve(Sig) x <- matrix(1, nrow=nrow(dados$coords)) return(drop(solve(t(x)%*%w%*%x, t(x)%*%w%*%dados$data))) } #parT<-(mu, s2, tau2, phi ) parTa <- c(22, 5, 13, 76) parTc <- c(25, 115, 0, 4) parTcb <- c(27, 170, 157, 33) parTic <- c(29, 8, 68, 109) res.a<-optim(c(log(sqrt(c(13,5))), log(76)), nloglik, dados=a, dmin=min(dist(a$coords)), method="B") c(mufun(c(exp(res.a$par[1:2])^2,exp(res.a$par[3])), a), c(exp(res.a$par[1:2])^2,exp(res.a$par[3]))) res.a res.c<-optim(c(log(sqrt(c(0.1,115))), log(4)), nloglik, dados=c, dmin=min(dist(c$coords)), method="S",control=list(maxit=100)) c(mufun(c(exp(res.c$par[1:2])^2,exp(res.c$par[3])), c), c(exp(res.c$par[1:2])^2,exp(res.c$par[3]))) res.c res.cb<-optim(c(log(sqrt(c(157,170))), log(33)), nloglik, dados=cb, dmin=min(dist(cb$coords)), method="S",control=list(maxit=100)) c(mufun(c(exp(res.cb$par[1:2])^2,exp(res.cb$par[3])), cb), c(exp(res.cb$par[1:2])^2,exp(res.cb$par[3]))) res.cb res.c<-optim(c(log(sqrt(c(13,5))), log(76)), nloglik, dados=c, dmin=min(dist(c$coords)), method="B") c(mufun(exp(res.c$par)^2, c),exp(res.c$par)^2) res.c nllik4.ind <- function(par, listadados, print.pars=TRUE) nloglik(par[1:4], listadados[[1]], print.pars) + nloglik(par[5:8], listadados[[2]], print.pars) + nloglik(par[9:12], listadados[[3]], print.pars) + nloglik(par[13:16], listadados[[4]], print.pars) nllik4.2a2 <- function(par, listadados, print.pars=TRUE) nloglik(par[1:4], listadados[[1]], print.pars) + nloglik(par[1:4], listadados[[2]], print.pars) + nloglik(par[5:8], listadados[[3]], print.pars) + nloglik(par[5:8], listadados[[4]], print.pars) nllik4.todos <- function(par, listadados, print.pars=TRUE) nloglik(par[1:4], listadados[[1]], print.pars) + nloglik(par[1:4], listadados[[2]], print.pars) + nloglik(par[1:4], listadados[[3]], print.pars) + nloglik(par[1:4], listadados[[4]], print.pars) #parT<-(mu, s2, tau2, phi ) parTa <- c(22, 5, 13, 76) parTc <- c(25, 115, 0, 4) parTcb <- c(27, 170, 157, 33) parTic <- c(29, 8, 68, 109) #par<-(tau, s ,phi) res.a<-optim(c(log(sqrt(c(3,5))), log(4)), nloglik, dados=a, dmin=min(dist(a$coords)), method="B") res.a <- optim(c(log(sqrt(c(13,5,76)))), nloglik, dados=a, dmin=min(dist(a$coords)), method="B") c(mufun(exp(res.a$par)^2, a),exp(res.a$par)^2) res.a res.c <- optim(log(sqrt(c(.1,100,4))), nloglik, dados=c, u=as.matrix(dist(c$coords)), dmin=min(dist(c$coords)), method="S",control=list(maxit=100)) c(mufun(exp(res.c$par)^2, c), exp(res.c$par)^2) res.c res.cb <- optim(log(sqrt(c(100,100,25))), nloglik, dados=cb, u=as.matrix(dist(cb$coords)), dmin=min(dist(cb$coords)), method="S",control=list(maxit=100)) c(mufun(exp(res.cb$par)^2, cb), exp(res.cb$par)^2) res.cb res.ic <- optim(log(sqrt(c(60,8,100))), nloglik, dados=ic, u=as.matrix(dist(ic$coords)), dmin=min(dist(ic$coords)), method="S",control=list(maxit=100)) c(mufun(exp(res.ic$par)^2, ic), exp(res.ic$par)^2) res.ic optim(c(ma, sa, ta, phia, mc, sc, tc, phic, mcb, scb, tcb, phicb, mic, sic, tic, phiic), nllik4.ind, listadados=list(a,c,cb,ic)) optim(c(maic,saic,taic,phiaic,mccb,sccb,tccb,phiccb), nllik4.2a2, listadados=list(a,ic, c,cb)) optim(c(mt,st,tt,phit), nllik4.todos, listadados=list(a,c,cb,ic)) #likfit(a, ini=c(sa, phia)) optim(c(mc, sc, tc, phic), nloglik, dados=c, method="L", lower=c(-Inf, -Inf, 0, 0)) # likfit(c, ini=c(sc, phic)) optim(c(mcb, scb, tcb, phicb), nloglik, dados=cb, method="L", lower=c(-Inf, -Inf, 0, 0)) optim(c(mic, sic, tic, phiic), nloglik, dados=ic, method="L", lower=c(-Inf, -Inf, 0, 0))