############################################################### # Normal Bivariada : # # Modelo Matern Multivariado Parcimonioso # ############################################################### 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 ## ############################## #d<-dist(as.matrix(gs)) #d #d<-as.vector(d) #d #length(d) 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) # ################################## ll<-function(pars, dados, logpars = FALSE, printpars = 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) if(conc){ if(logpars) pars[3:5] <- exp(pars[3:5]) if(is.null(names(pars))) names(pars) <- c("nu1","nu2","sig2","tau1","tau2","a","rho") parsL <- as.list(pars) if(printpars) print(unlist(parsL)) parsL$sig1 <- 1 } else{ if(logpars) pars[5:8] <- exp(pars[5:8]) if(is.null(names(pars))) names(pars) <- c("mu1","mu2", "nu1","nu2", "sig1","sig2","tau1","tau2","a","rho") parsL <- as.list(pars) if(printpars) print(unlist(parsL)) } if(evalq(abs(rho) >= sqrt(nu1 * nu2)/((nu1+nu2)/2), env=parsL, enclos=NULL)) return(-.Machine$double.xmax^0.5) SIGMA <- montaSigma(parsL, 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 <- 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 # ll<- -0.5*(dados$n*detSIGMA-HQ) ## aqui estava errado! 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) if(printpars) print(llik) #print(c(detSIGMA, HQ, llik)) return(llik) } ######################################### lista<-list(n=length(sim), n1=length(sim)/2, n2=length(sim)/2, X=sim, d=as.vector(dist(gs))) ## testando com varias formas de chamar a fun????o #function(pars, dados, logpars=FALSE, printpars=FALSE, conc = FALSE) ini<-c(mu1=1,mu2=2,nu1=1,nu2=1,sig1=sqrt(2),sig2=sqrt(2),tau1=sqrt(0.6),tau2=sqrt(0.5),a=1,rho=0.6) #debug(ll) ll(ini, dados=lista) #undebug(ll) ## com a reparametriza????o log() -- note que o valor da vero tem que ser o mesmo! ini<-c(mu1=1,mu2=2,nu1=1,nu2=1,sig1=log(sqrt(2)),sig2=log(sqrt(2)), tau1=log(sqrt(0.6)),tau2=log(sqrt(0.5)),a=1,rho=0.6) ll(ini, dados=lista,logpars=T) ## verossimilhanca concentrada ## ACHAR OS MU ini<-c(nu1=5,nu2=5,sig2=sqrt(2)/sqrt(2),tau1=sqrt(0.6)/sqrt(2),tau2=sqrt(0.5)/sqrt(2),a=1,rho=0.6) ll(ini, dados=lista, conc=TRUE) ## ini<-c(mu1=5,mu2=5,nu1=0.5,nu2=0.5,sig1=sqrt(2),sig2=sqrt(2),tau1=sqrt(0.5),tau2=sqrt(0.2),a=1/5,rho=0.5) ajuste <- optim(par=unname(ini), fn=ll, dados=lista, method="L-BFGS-B", lower=c(-Inf, -Inf, rep(1e-10, 7), -0.99), upper= c(rep(Inf, 9), 0.99), control=list(fnscale=-1),printpars=T) ## printpars=T, logpars=T ajuste[1:2] ###################################################################################### ## 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(mu1=5,mu2=5,nu1=1,nu2=1,sig1=sqrt(2),sig2=sqrt(2)/sqrt(2),tau1=sqrt(0.6)/sqrt(2),tau2=sqrt(0.5)/sqrt(2),a=1,rho=0.6) ajuste.inteiro <- optim(par=unname(ini), fn=ll, dados=lista, method="L-BFGS-B", lower=c(-Inf, -Inf, rep(1e-10, 7), -0.99), upper= c(rep(Inf, 9), 0.99), control=list(fnscale=-1),printpars=T) ajuste.inteiro 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(par=ini, 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 ajuste ##Utilizando a Verossimilhan??a concentrada: calculando mu1, mu2 e sigma1 pars.est<-list(nu1<-ajuste[[1]][[1]], nu2<-ajuste[[1]][[2]], sig2<-ajuste[[1]][[3]], tau1<-ajuste[[1]][[4]], tau2<-ajuste[[1]][[5]], a<-ajuste[[1]][[6]], rho<-ajuste[[1]][[7]], sig1<-1 ) names(pars.est) <- c('nu1','nu2','sig2','tau1','tau2','a','rho','sig1') dados<-list(n1=length(sim)/2, n2=length(sim)/2, X=sim, d=as.vector(dist(gs))) param<-function(pars.est,dados){ SIGMA <- montaSigma(pars.est, dados$d) tchSIGMA <- t(chol(SIGMA)) ONE <- kronecker(diag(1, 2), rep(1, dados$n1)) ONEchSIGMA <- solve(tchSIGMA, ONE) mus <- solve(crossprod(ONEchSIGMA), crossprod(ONEchSIGMA, solve(tchSIGMA, dados$X))) dif<-dados$X-ONE%*%mus sig1<-(1/(dados$n1+dados$n2))*(t(dif)%*%solve(tchSIGMA,dif)) sig2<-pars.est$sig2/sig1 tau1<-pars.est$tau1/sig1 tau2<-pars.est$tau2/sig1 nu1<-pars.est$nu1 nu2<-pars.est$nu2 a<-pars.est$a rho<-pars.est$rho parametros <- list(mus,sig1,sig2,tau1,tau2,nu1,nu2,a,rho) names(parametros) = c('mu','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') return(parametros) } mu.sig1<-param(pars.est,dados) mu.sig1