############################################################### # Normal Bivariada : # # Modelo Matern Multivariado Parcimonioso # # geral # ############################################################### rm(list=ls(all=TRUE)) #Bivariado: p=2, d=2 require(geoR) ########## ## Grid ## ########## set.seed(30) ## grid regular #gs <- expand.grid(seq(0,25, l=11), seq(0,25, l=11)) #gs ## grid irregular gs <- cbind(runif(121, 0, 25), runif(121, 0, 25)) gs plot(gs) ############################## ## distancia dos pontos ## ############################## dM <- dist(as.matrix(gs), upper=T, diag=T) dM #################### ## Simulando Dados## #################### ## valores fixados dos parametros para simulação LP <- list(nu1 = 0.7 , nu2 = 1.5, a = 1/4 , rho = 0.7 , sig1 = sqrt(0.9), sig2 = sqrt(0.5) , tau1 = sqrt(0.1), tau2 = sqrt(0.1) ) ## Duas opcoes para montar a matriz de covariâncias: ## 1. a partir do vetor com as distancias dos pares (triangulo inferior) ## 2. a partir da matrix com as distancias dos pares ## 1. #n=dimensao matriz lower2matrix <- function(x){ n <- (1+sqrt(1+8*length(x)))/2 M <- diag(n) M[lower.tri(M)] <- x M <- t(M) M[lower.tri(M)] <- x M } montaSigma <- function(parList, distPares){ evalq({ n <- (1+sqrt(1+8*length(distPares)))/2 ## condição/restrição if(abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2)) stop("restrições nos parâmetros não satisfeitas") c11 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=nu1) c11 <- sig1^2 * lower2matrix(c11) diag(c11) <- diag(c11) + tau1^2 c22 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=nu2) c22 <- sig2^2 * lower2matrix(c22) diag(c22) <- diag(c22) + tau2^2 c12 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=(nu1+nu2)/2) c12 <- rho * sig1 * sig2 * lower2matrix(c12) S <- rbind(cbind(c11, c12), cbind(c12, c22)) }, env=parList) } mS <- function(parList, distPares){ evalq({ ## condição/restrição if(abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2)) stop("restrições nos parâmetros não satisfeitas") c11 <- varcov.spatial(dists.l=distPares,cov.model="matern",cov.pars=c(sig1^2, 1/a), kappa=nu1, nugget=tau1^2)$varcov c22 <- varcov.spatial(dists.l=distPares,cov.model="matern",cov.pars=c(sig2^2, 1/a), kappa=nu2, nugget=tau2^2)$varcov c12 <- varcov.spatial(dists.l=distPares,cov.model="matern",cov.pars=c(rho*sig1*sig2, 1/a), kappa=(nu1+nu2)/2)$varcov S <- rbind(cbind(c11, c12), cbind(c12, c22)) }, env=parList) } ## 2. mSigma <- function(parList, distM){ local({ ## condição/restrição if(abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2)) stop("restrições nos parâmetros não satisfeitas") c11 <- cov.spatial(distM,cov.model="matern",cov.pars=c(sig1^2, 1/a),kappa=nu1) diag(c11) <- diag(c11) + tau1^2 c22 <- cov.spatial(distM,cov.model="matern",cov.pars=c(sig2^2, 1/a),kappa=nu2) diag(c22) <- diag(c22) + tau2^2 c12 <- cov.spatial(distM,cov.model="matern",cov.pars=c(rho * sig1 * sig2, 1/a),kappa=(nu1+nu2)/2) S <- rbind(cbind(c11, c12), cbind(c12, c22)) }, env=parList) } dV <- as.vector(dist(gs)) system.time(SIGMA <- montaSigma(LP, distPares=dV)) system.time(SIGMA1 <- mS(LP, distPares=dV)) dM <- as.matrix(dist(gs, diag=T,up=T)); dimnames(dM) <- list(NULL, NULL) system.time(SIGMA2 <- mSigma(LP, distM=dM)) all.equal(SIGMA, SIGMA1) all.equal(SIGMA, SIGMA2) ## o seguinte tb funciona: #names(LP) #LP$distPares <- as.vector(dist(gs)) #SIGMA <- montaSigma(LP) cholSIGMA <- chol(SIGMA) sim <- drop(crossprod(cholSIGMA, rnorm(dim(cholSIGMA)[1]))) sim length(sim) ## fazendo médias diferentes sim <- sim + c(rep(2, 121), rep(10, 121)) plot(variog(coords=gs, data=sim[1:121]), ty="b") lines(variog(coords=gs, data=sim[122:242]), ty="b", col=2) rm(SIGMA, SIGMA1, SIGMA2, cholSIGMA) ################################## # log(verrosimilhança) # ################################## geral<-function(pars1,pars2, printpars = FALSE, ...){ pars<-as.list(c(pars1,pars2)) rodaPARS <<- as.vector(pars) TODOS <- c("mu1","mu2", "nu1","nu2", "sig1","sig2","tau1","tau2","a","rho") if(any(is.na(match(TODOS, names(pars))))) stop("tá faltando parametro!!!") if(pars$nu1 <=0 || pars$nu2 <=0) return(-.Machine$double.xmax^0.5) # if(pars$sig1 <=0 || pars$sig2 <=0) return(-.Machine$double.xmax^0.5) ll<-function(pars, dados, logpars = FALSE, conc = FALSE) { # parsList recebe desvios (e não variâncias!) padrão com opção para log(desvio padrao) # dados : list(d, X, n1, n2, n)1e-10 if(!all(c("mu1","mu2","sig1") %in% names(pars1))) conc <- FALSE if(conc){ if(logpars) pars$sig2 <- exp(pars$sig2) pars$tau1 <- exp(pars$tau1) pars$tau2 <- exp(pars$tau2) pars$sig1 <- 1 } else{ if(logpars) pars$sig1 <- exp(pars$sig1) pars$sig2 <- exp(pars$sig2) pars$tau1 <- exp(pars$tau1) pars$tau2 <- exp(pars$tau2) } if(with(pars, abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2))) return(-.Machine$double.xmax^0.5) SIGMA <-montaSigma(pars, dados$d) tchSIGMA <- t(chol(SIGMA)) if(conc){ ## matrix (supondo n1 = n2) !!!!!! ONE <- kronecker(diag(1, 2), rep(1, dados$n1)) ## se n1 != n2: ## ONE <- cbind(c(rep(1, dados$n1), rep(0, dados$n2)), c(rep(0, dados$n1), rep(1, dados$n2))) ONEchSIGMA <- forwardsolve(tchSIGMA, ONE) mus <- solve(crossprod(ONEchSIGMA), crossprod(ONEchSIGMA, forwardsolve(tchSIGMA, dados$X))) } else mus <- c(pars$mu1,pars$mu2) #pars[1:2] Mu <- c(rep(mus[1],dados$n1), rep(mus[2],dados$n2)) res <- dados$X - Mu # HQ<-sum(res*drop(solve(SIGMA,res))) ## vamos fazer isto usando a dec. choleski z <- forwardsolve(tchSIGMA,res) HQ <- sum(z*z) detSIGMA <- 2*sum(log(diag(tchSIGMA))) ## det. calculado a partir da choleski if(conc) llik <- -0.5 * (dados$n * (log(2*pi) + log(HQ/dados$n)) + detSIGMA + dados$n) else llik<- -0.5 * (dados$n * log(2*pi) + detSIGMA + HQ) return(llik) } llik <- ll(pars, ...) if(printpars){ print(unlist(pars1)) print(llik) } return(llik) } ######################################### lista<-list(n=length(sim), n1=length(sim)/2, n2=length(sim)/2, X=sim, d=as.vector(dist(gs))) ini <-c(tau1=sqrt(0.6),tau2=sqrt(0.5),rho=0.6,nu1=1,nu2=1,sig2=sqrt(2),a=1) fixos <- c(mu1=1,mu2=2,sig1=sqrt(2)) ini <-c(c(mu1=1,mu2=2), log(c(sig1=sqrt(0.6),tau1=0.1,tau2=sqrt(0.5),sig2=sqrt(2))), rho=0.6,nu1=1,nu2=1,a=1) fixos=NULL #undebug(geral) geral(ini,pars2=fixos,dados=lista, conc=T) #undebug(geral) res1 <- optim(ini, geral, pars2=fixos, dados=lista,method="BFGS", control=list(fnscale=-1, maxit=6000), conc=T, printpars=T, logp=T) res1 res2<-optim(ini, geral, pars2=fixos, dados=lista,method="L-BFGS-B", lower=c(-Inf,-Inf,1e-5,1e-5,1e-5,1e-5, -1,rep(0, 3)), upper= c(rep(Inf,6),1,rep(Inf, 3)), control=list(fnscale=-1, maxit=3000), conc=T, printpars=T) res2 ## ini<-c(mu1=1,mu2=2,nu1=1,nu2=1,sig1=sqrt(2),sig2=sqrt(2),a=1,tau1=sqrt(0.6),tau2=sqrt(0.5),rho=0.6) ajuste <- optim(ini, fn=geral,pars2=NULL, dados=lista, method="L-BFGS-B", lower=c(-Inf, -Inf, rep(0, 7), -1), upper= c(rep(Inf, 9), 1), control=list(fnscale=-1,maxit=3000)) ## testar com diferentes valores iniciais e ver se converge para os mesmos lugares ## verificar tb opcoes com e sem log ## note que com log ou reparametrizacoes nao precisa de L-BFGS-B !! ## valores iniciais devem ser mais ou menos compativeis com: ## mu1: media da primeira variavel ## mu2: media da segunda variavel ## tau1^2 + sig1^2 : variancia da primeira variavel ## tau2^2 + sig2^2 : variancia da segunda variavel ini<-c(nu1=1,nu2=1,sig2=sqrt(2)/sqrt(2),tau1=sqrt(0.6)/sqrt(2),tau2=sqrt(0.5)/sqrt(2),a=1,rho=0.6) ajuste <- optim(pars1=ini,pars2=fixo ,fn=ll, dados=lista, method="L-BFGS-B", lower=c(rep(0, 6), -1), upper= c(rep(Inf, 6), 1), control=list(fnscale=-1), conc=T, printpars=T) # logpars=T