# ESALQ/USP # GEOESTATÍSTICA, 2012 # EXERCÍCIO: Estimação de parâmetros do semivariograma utilizando # a distribuição t-multivariada # ANDERSON RODRIGO DA SILVA # Utilizando os dados 'ca20' da geoR require(geoR) data(ca20) v <- variog(ca20, max.dist=700) plot(v) ## ----------------------------------------------------------------------- ## Estimação por ML, considerando a distribuição t multivariada ## ----------------------------------------------------------------------- # install.packages("mvtnorm") require(mvtnorm) ltheta <- function(mu, tausq, sigmasq, phi) { n <- nrow(ca20$coords); I <- diag(n) H <- as.matrix(dist(ca20$coords, upper=TRUE, diag=TRUE)) S <- tausq * I + sigmasq * exp(-H/phi) return(-dmvt(x=ca20$data, delta=rep(mu, n), sigma=S, df = 1, log = TRUE)) } ltheta(40, 20, 100, 100) t.est <- mle(minuslogl=ltheta, start=list(mu=40, tausq=20, sigmasq=100, phi=100), method="L-BFGS-B", lower=c(0,0,0,0)) t.est ## ----------------------------------------------------------------------- ## Estimação por ML, considerando a distribuição Gaussiana ## ----------------------------------------------------------------------- # install.packages("mnormt") require(mnormt) ltheta2 <- function(mu, tausq, sigmasq, phi) { n <- nrow(ca20$coords); I <- diag(n) H <- as.matrix(dist(ca20$coords, upper=TRUE, diag=TRUE)) S <- tausq * I + sigmasq * exp(-H/phi) return(-dmnorm(x=ca20$data, mean=rep(mu, n), varcov=S, log=TRUE)) } ltheta2(40,20,100,100) norm.est <- mle(minuslogl=ltheta2, start=list(mu=40, tausq=20, sigmasq=100, phi=100), method="L-BFGS-B", lower=c(0,0,0,0)) norm.est lik <- likfit(ca20, ini=c(100,100)) # com a geoR lik # Comparando resultados summary(t.est) summary(norm.est) # observa-se que as variâncias das estimativas dos parâmetros de escala (tausq e sigmasq) # são maiores quando se utiliza a distribuição t-multivariada. # verificando os ajustes ao semivariograma experimental plot(v) lines(lik) d <- seq(1:900) points(t.est@coef[2] + t.est@coef[3]*(1-exp(-d/t.est@coef[4])) ~ d, type="l", col="red") legend(x=0, y=150, c("t-Student", "Gaussiana"), col=c(2,1), lty=c(1,1), cex=0.8) #