############################################################### # Normal Bivariada : # # Modelo Matern Multivariado Parcimonioso # ############################################################### rm(list=ls(all=TRUE)) #Bivariado: p=2, d=2 require(geoR) require(fields) require(RandomFields) require(MASS) ########## ## Grid ## ########## set.seed(29) gs <- expand.grid(seq(0,30, l=11), seq(0,30, l=11)) gs ############################## ## distancia dos pontos ## ############################## d<-dist(as.matrix(gs)) d d<-as.vector(d) d length(d) #################### ## Simulando Dados## #################### #n=dimensao matriz lower2matrix <- function(x){ n <- (1+sqrt(1+8*length(x)))/2 M <- diag(n) M[lower.tri(M)] <- x M <- M + t(M) - diag(n) } nu1 <- 1 nu2 <- 1 a <- 1/4 phi<-4 rho <- 0.7 sig1<-0.6 sig2<-0.5 #matriz de covariância(SIGMA) cov11<-cov.spatial(d,cov.model="matern",cov.pars=c(1,phi),kappa=nu1) c11 <- sig1 * lower2matrix(cov11)+0.1*diag(length(gs$Var1)) cov22<-cov.spatial(d,cov.model="matern",cov.pars=c(1,phi),kappa=nu2) c22<-sig2*lower2matrix(cov22)+0.1*diag(length(gs$Var2)) cov12<-cov.spatial(d,cov.model="matern",cov.pars=c(1,phi),kappa=(nu1+nu2)/2) c12<-rho*sqrt(sig1)*sqrt(sig2)*lower2matrix(cov12) ## condição abs(rho) < sqrt(nu1 * nu2)/((nu1+nu2)/2) SIGMA<-rbind(cbind(c11, c12), cbind(t(c12), c22)) cholSIGMA<-chol(SIGMA) sim <- as.vector(crossprod(cholSIGMA) %*% rnorm(dim(SIGMA)[1])) sim ################################## # 1 # ################################## ll<-function(pars,dados) { #pars=c(mu1,mu2,nu1,nu2,sig2.1,sig2.2,tau1,tau2,a,rho12),dados=list(n,n1,n2,y,d) cov11<-cov.spatial(dados$d,cov.model="matern",cov.pars=c(1,pars[9]),kappa=pars[3]) c11 <- pars[5] * lower2matrix(cov11)+pars[7]*diag(length(dados$n1)) cov22<-cov.spatial(d,cov.model="matern",cov.pars=c(1,pars[9]),kappa=pars[4]) c22<-pars[6]*lower2matrix(cov22)+pars[8]*diag(length(dados$n2)) cov12<-cov.spatial(d,cov.model="matern",cov.pars=c(1,pars[9]),kappa=(pars[3]+pars[4])/2) c12<-pars[12]*sqrt(pars[5])*sqrt(pars[6])*lower2matrix(cov12) SIGMA<-rbind(cbind(c11, c12), cbind(t(c12), c22)) detSIGMA<-determinant(SIGMA,log=T)$modulus #JA APLICANDO LOG Mu<-c(rep(pars[1],dados[2]),rep(pars[2],dados[3])) Mu<-as.matrix(Mu,ncol=1) res<-dados$y - Mu HQ<-sum(res*solve(SIGMA,res)) ll<- -0.5*(dados$n*detSIGMA-HQ) return(ll) } ################################## # 2 # ################################## ll<-function(pars,dados) { names(pars) <- c("mu1","mu2","nu1","nu2","sig2.1","sig2.2","tau1","tau2","a","rho12") print(pars) a=as.numeric(pars[9]);nu1=as.numeric(pars[3]);sig2.1=as.numeric(pars[5]);tau1=as.numeric(pars[7]); nu2=as.numeric(pars[4]);sig2.2=as.numeric(pars[6]);tau2=as.numeric(pars[8]); mu1=as.numeric(pars[1]);mu2=as.numeric(pars[2]);rho12=as.numeric(pars[10]) if(abs(rho) >= sqrt(nu1 * nu2)/((nu1+nu2)/2)) return(-.Machine$double.xmax^0.5) #if(abs(rho) >= sqrt(nu1 * nu2)/((nu1+nu2)/2)) return(NA) cov11<-cov.spatial(dados$d,cov.model="matern",cov.pars=c(1,a),kappa=nu1) c11 <- sig2.1 * lower2matrix(cov11)+tau1*diag(dados$n1) cov22<-cov.spatial(dados$d,cov.model="matern",cov.pars=c(1,a),kappa=nu2) c22<-sig2.2*lower2matrix(cov22)+tau2*diag(dados$n2) cov12<-cov.spatial(dados$d,cov.model="matern",cov.pars=c(1,a),kappa=(nu1+nu2)/2) c12<-rho12*sqrt(sig2.1)*sqrt(sig2.2)*lower2matrix(cov12) SIGMA<-rbind(cbind(c11, c12), cbind(t(c12), c22)) detSIGMA<-determinant(SIGMA,log=T)$modulus #JA APLICANDO LOG attributes(detSIGMA) <- NULL Mu<- c(rep(mu1,dados$n1),rep(mu2,dados$n2)) Mu<-as.matrix(Mu,ncol=1) res<-dados$X-Mu HQ<-sum(res*solve(SIGMA,res)) llik<- -0.5*(dados$n*detSIGMA-HQ) print(llik) return(llik) } ######################################### lista<-list(n=242, n1=121, n2=121, X=sim, d=as.vector(dist(gs))) ini<-c(mu1=1,mu2=2,nu1=1,nu2=1,sig2.1=2,sig2.2=2,tau1=0.6,tau2=0.5,a=0.8,rho12=0.6) ll(ini, dados=lista) ajuste <- optim(par=ini, fn=ll, dados=lista, method="L-BFGS-B", lower=c(0,0,0,0,0,0,0,0,0,-1), upper= c(rep(Inf, 9), 1), control=list(fnscale=-1)) ajuste[1]