################################################################################ # ESTUDO DE SIMULAÇÃO DE MODELOS GEOESTATÍSTICO BIVARIADOS # # Simulação, Estimação MV, Krigagem e análise de resultados # # # # Bruno Henrique Fernandes Fonseca e Paulo Justiniano Ribeiro Jr. # # Piracicaba, 2008 # ################################################################################ ############################# ANÁLISE DO BGCCM ################################# ## Função para simular dados com o BGCCM require(MASS) require(geoR) rBGCCM_sim <- function(n1,n2,N,n_arm,perc_co_loc=1,medias=c(0,0),sigmas,phis, kappas=c(.5,.5,.5),sem_grid=1,sem_cam_ale=1){ set.seed((sem_grid*3331)) grx1 <- runif(n1) set.seed((sem_grid*3771)) gry1 <- runif(n1) gr1 <- unname(cbind(grx1[1:n1],gry1[1:n1])) set.seed((sem_grid*3713)) grx2 <- runif(n2) set.seed((sem_grid*3317)) gry2 <- runif(n2) gr2 <- unname(cbind(grx2,gry2)) set.seed((sem_grid*3137)) grx3 <- runif(n_arm) set.seed((sem_grid*3337)) gry3 <- runif(n_arm) gr3 <- unname(cbind(grx3,gry3)) n <- n1+n2 n_co <- round(min(n1,n2)*perc_co_loc,0) k1 <- c((n1+n_co+1):n) k2 <- c(1:(n2-n_co)) k1_1 <- c((n_co+1):n1) k2_1 <- c(1:(n1-n_co)) k3 <- c(1:n) k4 <- c(1:n1) k4_1 <- c(1:(n1+n_arm)) k5 <- c((n1+1):n) k5_1 <- c((n1+n_arm+1):(n+2*n_arm)) k5_2 <- c((n1+n_arm+1):(n1+n_arm+n2)) k6 <- c(1:n2) k6_1 <- c(1:(n2+n_arm)) k7 <- c((n1+n_arm+n_co+1):(n+n_arm)) k8 <- c(1:(n+2*n_arm)) if(n1==n2) gr <- rbind(gr1,gr3,gr1,gr3) if(n1==n2 & perc_co_loc < 1) gr[k7,] <- gr2[k2,] if(n1>n2) gr <- rbind(gr1,gr3,gr1[k6,],gr3) if(n1>n2 & perc_co_loc < 1) gr[k7,] <- gr2[k2,] if(n2>n1) gr <- rbind(gr2[k4,],gr3,gr2,gr3) if(n2>n1 & perc_co_loc < 1) gr[k1_1,] <- gr1[k2_1,] m_1 <- medias[1] m_2 <- medias[2] s01 <- sigmas[1] s02 <- sigmas[3] s1 <- sigmas[2] s2 <- sigmas[4] phi0 <- phis[1] phi1 <- phis[2] phi2 <- phis[3] kappa0 <- kappas[1] kappa1 <- kappas[2] kappa2 <- kappas[3] set.seed(sem_cam_ale*31) R01 <- grf(n=(n1+n_arm),grid=gr[k4_1,],nsim=N,cov.pars=c(s01^2,phi0),kappa=kappa0) set.seed(sem_cam_ale*31) R02 <- grf(n=(n2+n_arm),grid=gr[k5_1,],nsim=N,cov.pars=c(s02^2,phi0),kappa=kappa0) set.seed(sem_cam_ale*37) R1 <- grf(n=(n1+n_arm),grid=gr[k4_1,],nsim=N,cov.pars=c(s1^2,phi1),kappa=kappa1) set.seed(sem_cam_ale*71) R2 <- grf(n=(n2+n_arm),grid=gr[k5_1,],nsim=N,cov.pars=c(s2^2,phi2),kappa=kappa2) y1 <- matrix(0,(n1+n_arm),N) y2 <- matrix(0,(n2+n_arm),N) if(N >1){ y1 <- m_1 + R01$data + R1$data y2 <- m_2 + R02$data + R2$data } if(N == 1){ y1[,1] <- m_1 + R01$data + R1$data y2[,1] <- m_2 + R02$data + R2$data } saida <- list() geo1 <- cbind(gr[k4,],y1[k4,]) geo2 <- cbind(gr[k5_2,],y2[k6,]) geo_arm1 <- cbind(gr3,y1[(n1+1):(n1+n_arm),]) geo_arm2 <- cbind(gr3,y2[(n2+1):(n2+n_arm),]) saida$geo_arm1 <- as.geodata(geo_arm1,data.col=c(3:(2+N))) saida$geo_arm2 <- as.geodata(geo_arm2,data.col=c(3:(2+N))) saida$geo1 <- as.geodata(geo1,data.col=c(3:(2+N))) saida$geo2 <- as.geodata(geo2,data.col=c(3:(2+N))) saida$geo1$coords <- unname(saida$geo1$coords) saida$geo2$coords <- unname(saida$geo2$coords) saida$geo_arm1$coords <- unname(saida$geo_arm1$coords) saida$geo_arm2$coords <- unname(saida$geo_arm2$coords) saida$teste1 <- R01$data saida$teste2 <- R02$data return(saida) } ## Função para fazer estimação por MV dos parâmetros do BGCCM emvBGCCM_sim <- function(geo1,geo2,corr_ini,ruido=F,kappas=c(0.5,0.5,0.5),metodo="Nelder-Mead",inferior=-Inf,superior=Inf, controle=list(),hessiana=F){ n1 <- nrow(geo1$coords) n2 <- nrow(geo2$coords) n <- n1+n2 N <- ncol(geo1$data) if(is.null(N)) N <- 1 k1 <- c(1:n1) k2 <- c((n1+1):n) k3 <- c(1:n) k4 <- c(1:n2) dados <- matrix(0,n,N) if(N > 1) for(i in 1:N) dados[,i] <- c(geo1$data[,i],geo2$data[,i]) if(N == 1) dados[,1] <- c(geo1$data,geo2$data) dists <- unname(as.matrix(dist(rbind(geo1$coords,geo2$coords),diag=TRUE,upper=TRUE))) X <- matrix(c(rep(1,n1),rep(0,n2),rep(0,n1),rep(1,n2)),n,2) kappa0 <- kappas[1] kappa1 <- kappas[2] kappa2 <- kappas[3] ctes <- list() ctes$n1 <- n1 ctes$n2 <- n2 ctes$dados <- dados ctes$dists <- dists ctes$X <- X est_num<-function(p,Ct,rui=F,kap,repet){ if ( any(p < 1e-20)) return(.Machine$double.xmax^.5) if ( any(p[4:6] > 1)) return(.Machine$double.xmax^.5) if ( any(p > 100)) return(.Machine$double.xmax^.5) R0 <- (((2^(kap[1]-1))*gamma(kap[1]))^(-1)) * ((Ct$dists/p[4])^(kap[1])) * besselK(Ct$dists/p[4],nu=kap[1]) R0[ Ct$dists < 0.005 ] <- 1 R1 <- (((2^(kap[2]-1))*gamma(kap[2]))^(-1)) * ((Ct$dists/p[5])^(kap[2])) * besselK(Ct$dists/p[5],nu=kap[2]) R1[ Ct$dists < 0.005 ] <- 1 R2 <- (((2^(kap[3]-1))*gamma(kap[3]))^(-1)) * ((Ct$dists/p[6])^(kap[3])) * besselK(Ct$dists/p[6],nu=kap[3]) R2[ Ct$dists < 0.005 ] <- 1 V1 <- R0[1:Ct$n1,1:Ct$n1] + p[2]^2 * R1[1:Ct$n1,1:Ct$n1] V2 <- p[1]^2 * R0[(Ct$n1+1):(Ct$n1+Ct$n2),(Ct$n1+1):(Ct$n1+Ct$n2)] + p[3]^2 * R2[(Ct$n1+1):(Ct$n1+Ct$n2),(Ct$n1+1):(Ct$n1+Ct$n2)] V3 <- p[1] * R0 [1:Ct$n1,(Ct$n1+1):(Ct$n1+Ct$n2)] V <- cbind(rbind(V1,t(V3)),rbind(V3,V2)) if(rui == TRUE) {diag(V) <- diag(V) + c(rep((p[7]^2),Ct$n1),rep((p[8]^2),Ct$n2))} Vinv <- ginv(V) tVinvX <- crossprod(Vinv,Ct$X) Q_hat <- drop(sum(crossprod(Ct$dados[,repet],Vinv)*Ct$dados[,repet]) - crossprod(Ct$dados[,repet],tVinvX)%*%ginv(crossprod(tVinvX,Ct$X))%*%crossprod(tVinvX,Ct$dados[,repet])) lQ_hat <- log(Q_hat) ldv <- log(det(V)) if(is.infinite(lQ_hat) | is.nan(lQ_hat) | is.infinite(ldv) | is.nan(ldv)) return(.Machine$double.xmax^.5) ll <- drop(-0.5*((Ct$n1+Ct$n2)*(log(2*pi)+lQ_hat-log(Ct$n1+Ct$n2)+1)+ldv)) if(is.infinite(ll) | is.nan(ll)) return(.Machine$double.xmax^.5) return(-ll) } resultados <- list() resultados$betas <- matrix(0,N,2) colnames(resultados$betas) <- c("mu1","mu2") resultados$sigmas <- matrix(0,N,4) colnames(resultados$sigmas) <- c("s01","s1","s02","s2") resultados$phis <- matrix(0,N,3) colnames(resultados$phis) <- c("phi0","phi1","phi2") if(ruido == TRUE){ resultados$taus <- matrix(0,N,2) colnames(resultados$taus) <- c("tau1","tau2") resultados$re_parametros <- matrix(0,N,5) colnames(resultados$re_parametros) <- c("eta","nu1","nu2","psi1","psi2") if(hessiana == TRUE){ resultados$hessiana <- matrix(0,(8*N),8) colnames(resultados$hessiana) <- c("eta","nu1","nu2","phi0","phi1","phi2","psi1","psi2") rownames(resultados$hessiana) <- rep(c("eta","nu1","nu2","phi0","phi1","phi2","psi1","psi2"),N) } } if(ruido==FALSE){ resultados$re_parametros <- matrix(0,N,3) colnames(resultados$re_parametros) <- c("eta","nu1","nu2") if(hessiana == TRUE){ resultados$hessiana <- matrix(0,(6*N),6) colnames(resultados$hessiana) <- c("eta","nu1","nu2","phi0","phi1","phi2") rownames(resultados$hessiana) <- rep(c("eta","nu1","nu2","phi0","phi1","phi2"),N) } } resultados$log_lik <- matrix(0,N,1) resultados$kappas <- c(kappa0=kappas[1],kappa1=kappas[2],kappa2=kappas[3]) resultados$erros <- matrix(0,N,7) colnames(resultados$erros) <- c("repet","eta","nu1","nu2","phi0","phi1","phi2") convergencia <- rep(0,N) for(i in 1:N){ var1_amost <- var(geo1$data[,i]) var2_amost <- var(geo2$data[,i]) var_ini <- NULL var_ini[1] <- sqrt(0.8*var1_amost) var_ini[2] <- sqrt(0.2*var1_amost) var_ini[3] <- sqrt(0.8*var2_amost) var_ini[4] <- sqrt(0.2*var2_amost) repar_ini <- 0 repar_ini[1] <- var_ini[3]/var_ini[1] repar_ini[2] <- var_ini[2]/var_ini[1] repar_ini[3] <- var_ini[4]/var_ini[1] repar_ini[4] <- corr_ini[1] repar_ini[5] <- corr_ini[2] repar_ini[6] <- corr_ini[3] if(ruido == TRUE) { repar_ini[7] <- var_ini[5]/var_ini[1] repar_ini[8] <- var_ini[6]/var_ini[1] } estimativas <- try(optim(p=repar_ini,fn=est_num,Ct=ctes,rui=ruido,kap=c(kappa0,kappa1,kappa2),repet=i, method=metodo,lower=inferior,upper=superior,control=controle,hessian=hessiana)) if(class(estimativas) == "try-error"){ resultados$erros[i,] <- c(i,NaN,NaN,NaN,NaN,NaN,NaN) resultados$betas[i,] <- c(NaN,NaN) resultados$sigmas[i,] <- c(NaN,NaN,NaN,NaN) resultados$phis[i,] <- c(NaN,NaN,NaN) if(ruido == TRUE){ resultados$taus[i,] <- c(NaN,NaN) resultados$re_param[i,] <- c(NaN,NaN,NaN,NaN,NaN) } if(ruido==FALSE) resultados$re_parametros[i,] <- c(NaN,NaN,NaN) resultados$log_lik[i,] <- NaN if(hessiana == TRUE){ resultados$hessiana[((length(estimativas$par)*i)-(length(estimativas$par)-1)):(i*length(estimativas$par)),] <- NaN } } if(class(estimativas) != "try-error"){ if(estimativas$convergence == 0){ convergencia[i] <- 1 } print(i) R0 <- (((2^(kappa0-1))*gamma(kappa0))^(-1)) * ((dists/estimativas$par[4])^(kappa0)) * besselK(dists/estimativas$par[4],nu=kappa0) R0[dists < 0.005] <- 1 R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/estimativas$par[5])^(kappa1)) * besselK(dists/estimativas$par[5],nu=kappa1) R1[dists < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists/estimativas$par[6])^(kappa2)) * besselK(dists/estimativas$par[6],nu=kappa2) R2[dists < 0.005] <- 1 V1 <- R0[k1,k1] + estimativas$par[2]^2 * R1[k1,k1] V2 <- estimativas$par[1]^2 * R0[k2,k2] + estimativas$par[3]^2 * R2[k2,k2] V3 <- estimativas$par[1] * R0[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) VinvX <- ginv(V)%*%X mu <- drop(ginv(crossprod(VinvX,X))%*%crossprod(VinvX,dados[,i])) dif_y_mu <- dados[,i]-(X%*%mu) sigmasq <- drop((1/n) * crossprod(dif_y_mu,ginv(V)) %*% dif_y_mu) #if(sigmasq < 0) sigmasq <- drop((1/n) * crossprod(dif_y_mu,chol2inv(V)) %*% dif_y_mu) if(sigmasq < 0) resultados$erros[i,] <- c(i,estimativas$par[1],estimativas$par[2],estimativas$par[3],estimativas$par[4],estimativas$par[5],estimativas$par[6]) sigma01 <- sqrt(sigmasq) sigma02 <- estimativas$par[1] *sigma01 sigma1 <- estimativas$par[2] * sigma01 sigma2 <- estimativas$par[3] * sigma01 if(ruido == TRUE) { tau_1 <- estimativas$par[7]*sigma01 tau_2 <- estimativas$par[8]*sigma01 } resultados$betas[i,] <- mu resultados$sigmas[i,] <- c(sigma01,sigma1,sigma02,sigma2) resultados$phis[i,] <- c(estimativas$par[4],estimativas$par[5],estimativas$par[6]) if(ruido == TRUE){ resultados$taus[i,] <- c(tau_1,tau_2) resultados$re_param[i,] <- c(estimativas$par[1],estimativas$par[2],estimativas$par[3],estimativas$par[7],estimativas$par[8]) } if(ruido==FALSE) resultados$re_parametros[i,] <- c(estimativas$par[1],estimativas$par[2],estimativas$par[3]) resultados$log_lik[i,] <- -estimativas$value if(hessiana == TRUE){ resultados$hessiana[((length(estimativas$par)*i)-(length(estimativas$par)-1)):(i*length(estimativas$par)),] <- estimativas$hessian } } } resultados$erros <- subset(resultados$erros,resultados$erros[,"repet"]>0) print(c("convergência",convergencia)) class(resultados) <- "emvBGCCM" return(resultados) } ## Função para fazer krigagem com o BGCCM predBGCCM_sim <- function(geos,parametros){ geo1 <- geos$geo1 geo2 <- geos$geo2 geo_a1 <- geos$geo_arm1 geo_a2 <- geos$geo_arm2 if(nrow(parametros$erros)>0){ geo1$data <- geo1$data[,-c(parametros$erros[,1])] geo2$data <- geo2$data[,-c(parametros$erros[,1])] geo_a1$data <- geo_a1$data[,-c(parametros$erros[,1])] geo_a2$data <- geo_a2$data[,-c(parametros$erros[,1])] } grp <- geo_a1$coords coords1 <- geo1$coords coords2 <- geo2$coords n1 <- nrow(coords1) n2 <- nrow(coords2) n <- n1+n2 npr <- nrow(grp) nt <- npr + n k1 <- c(1:n1) k2 <- c(1:n2) k3 <- c(1:npr) k4 <- c((npr+1):(npr+n1)) k5 <- c((npr+1):(2*npr)) k6 <- c((npr+n1+1):nt) k7 <- c((n1+1):n) k8 <- c((npr+1):nt) k9 <- c(((2*npr)+1):((2*npr)+n)) N <- ncol(geo1$data) if(is.null(N)) N <- 1 dados <- matrix(0,n,N) if(N > 1) for(i in 1:N) dados[,i] <- c(geo1$data[,i],geo2$data[,i]) if(N == 1) dados[,1] <- c(geo1$data,geo2$data) gr_dados <- rbind(coords1,coords2) gr <- rbind(grp,gr_dados) dists <- matrix(0,nt,n) for(i in 1:nt){ dists[i,] <- spDistsN1(gr_dados,gr[i,]) } X <- matrix(c(rep(1,npr),rep(0,npr),rep(1,n1),rep(0,n2),rep(0,npr),rep(1,npr),rep(0,n1),rep(1,n2)),(2*npr+n),2) betas <- parametros$betas sigmas <- parametros$sigmas phis <- parametros$phis if(nrow(parametros$erros)>0){ betas <- matrix(c(parametros$betas[-c(parametros$erros[,1]),]),N,2) sigmas <- matrix(c(parametros$sigmas[-c(parametros$erros[,1]),]),N,4) phis <- matrix(c(parametros$phis[-c(parametros$erros[,1]),]),N,3) } if(!is.null(parametros$taus)) taus <- parametros$taus kappas <- parametros$kappas saida <- list() saida$pred1 <- matrix(0,npr,N) saida$pred2 <- matrix(0,npr,N) saida$var_pred1 <- matrix(0,npr,N) saida$var_pred2 <- matrix(0,npr,N) for(l in 1:N){ R0 <- (((2^(kappas[1]-1))*gamma(kappas[1]))^(-1)) * ((dists/phis[l,1])^(kappas[1])) * besselK(dists/phis[l,1],nu=kappas[1]) R0[dists < 0.005] <- 1 R1 <- (((2^(kappas[2]-1))*gamma(kappas[2]))^(-1)) * ((dists[k4,k1]/phis[l,2])^(kappas[2])) * besselK(dists[k4,k1]/phis[l,2],nu=kappas[2]) distaux <- dists[k4,k1] R1[distaux < 0.005] <- 1 R2 <- (((2^(kappas[3]-1))*gamma(kappas[3]))^(-1)) * ((dists[k6,k7]/phis[l,3])^(kappas[3])) * besselK(dists[k6,k7]/phis[l,3],nu=kappas[3]) distaux <- dists[k6,k7] R2[distaux < 0.005] <- 1 Sig1 <- (sigmas[l,1]^2)*R0[k4,k1] + (sigmas[l,2]^2)*R1 Sig2 <- (sigmas[l,3]^2)*R0[k6,k7] + (sigmas[l,4]^2)*R2 Sig12 <- (sigmas[l,1]*sigmas[l,3])*R0[k4,k7] Sig <- cbind(rbind(Sig1,t(Sig12)),rbind(Sig12,Sig2)) if(!is.null(parametros$taus)){ diag(Sig) <- diag(Sig) + c(rep(taus[l,1],n1),rep(tau[l,2],n2)) } R1c <- (((2^(kappas[2]-1))*gamma(kappas[2]))^(-1)) * ((dists[k3,k1]/phis[l,2])^(kappas[2])) * besselK(dists[k3,k1]/phis[l,2],nu=kappas[2]) distaux <- dists[k3,k1] R1c[distaux < 0.005] <- 1 Sig_c1 <- (sigmas[l,1]^2)*R0[k3,k1] + (sigmas[l,2]^2)*R1c Sig_c12 <- (sigmas[l,1]*sigmas[l,3]) * R0[k3,k7] Sig_c <- cbind(Sig_c1,Sig_c12) m_p <- as.vector(X%*%betas[l,]) m_z_dadoy <- m_p[k3] - (Sig_c%*%ginv(Sig)%*%(m_p[k9]-dados[,l])) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred1[,l] <- m_z_dadoy saida$var_pred1[,l] <- ((sigmas[l,1]^2) + (sigmas[l,2]^2)) - diag(pondera_var) R2c <- (((2^(kappas[3]-1))*gamma(kappas[3]))^(-1)) * ((dists[k3,k7]/phis[l,3])^(kappas[3])) * besselK(dists[k3,k7]/phis[l,3],nu=kappas[3]) distaux <- dists[k3,k7] R2c[distaux < 0.005] <- 1 Sig_c2 <- (sigmas[l,3]^2)*R0[k3,k7] + (sigmas[l,4]^2)*R2c Sig_c21 <- (sigmas[l,1]*sigmas[l,3]) * R0[k3,k1] Sig_c <- cbind(Sig_c21,Sig_c2) m_z_dadoy <- m_p[k5] - (Sig_c%*%ginv(Sig)%*%(m_p[k9]-dados[,l])) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred2[,l] <- m_z_dadoy saida$var_pred2[,l] <- ((sigmas[l,3]^2) + (sigmas[l,4]^2)) - diag(pondera_var) } return(saida) } ## Simulando sim1 <- rBGCCM_sim(n1=100,n2=100,N=5,n_arm=20, perc_co_loc=1,medias=c(150,60),sigmas=c(8,4,5,2),phis=c(0.25,0.20,0.20), kappas=c(0.5,0.5,0.5),sem_cam_ale=1) ## Estimando est1 <- emvBGCCM_sim(sim1$geo1,sim1$geo2,corr_ini=c(0.25,0.2,0.2),kappas=c(0.5,0.5,0.5), hessiana=T,metodo="Nelder-Mead",controle=list(maxit=10000)) ## Krigando pr1 <- predBGCCM_sim(geos=sim1,parametros=est1) ## Análise # Resultados da estimação betas5 <- est1$betas sigmas5 <- est1$sigmas phis5 <- est1$phis re_parametros5 <- est1$re_parametros hessianas5 <- est1$hessiana log_liks5 <- est1$log_lik # Resultados da predição krigs1 <- pr1$pred1 krigs2 <- pr1$pred2 var_krigs1 <- pr1$var_pred1 var_krigs2 <- pr1$var_pred2 # Preparando constantes X <- matrix(c(rep(1,100),rep(0,100),rep(0,100),rep(1,100)),ncol=2) dados1 <- sim1$geo1$data dados2 <- sim1$geo2$data dados <- rbind(dados1,dados2) dists <- unname(as.matrix(dist(rbind(sim1$geo1$coords,sim1$geo2$coords),diag=TRUE,upper=TRUE))) kappa0 <- 0.5 kappa1 <- 0.5 kappa2 <- 0.5 k1 <- c(1:100) k2 <- c(101:200) # VARIÂNCIAS ESTIMADAS DOS BETAS e DO SIGMA = SIGMA01 var.betas <- matrix(0,nrow=nrow(betas5),ncol=2) var.s01 <- rep(0,nrow(sigmas5)) mat.beta.s01 <- matrix(0,nrow=3*nrow(betas5),3) for(i in 1:nrow(betas5)){ R0 <- (((2^(kappa0-1))*gamma(kappa0))^(-1)) * ((dists/phis5[i,1])^(kappa0)) * besselK(dists/phis5[i,1],nu=kappa0) R0[dists < 0.005] <- 1 R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/phis5[i,2])^(kappa1)) * besselK(dists/phis5[i,2],nu=kappa1) R1[dists < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists/phis5[i,3])^(kappa2)) * besselK(dists/phis5[i,3],nu=kappa2) R2[dists < 0.005] <- 1 V1 <- R0[k1,k1] + re_parametros5[i,2]^2 * R1[k1,k1] V2 <- re_parametros5[i,1]^2 * R0[k2,k2] + re_parametros5[i,3]^2 * R2[k2,k2] V3 <- re_parametros5[i,1] * R0[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) iV <- ginv(V) Q_hat <- t(dados[,i] - X%*%betas5[i,]) %*% iV %*% (dados[,i] - X%*%betas5[i,]) i.mat.beta.s01 <- matrix(0,3,3) i.mat.beta.s01[1:2,1:2] <- (t(X)%*%iV%*%X)/sigmas5[i,1]^2 i.mat.beta.s01[3,3] <- (2*nrow(dados)^2)/Q_hat i.mat.beta.s01[1:2,3] <- (2*(t(X)%*%iV%*%dados[,i] - t(X)%*%iV%*%X%*%betas5[i,])) / sigmas5[i,1]^3 i.mat.beta.s01[3,1:2] <- t(i.mat.beta.s01[1:2,3]) mat.beta.s01[(3*i-2):(3*i),] <- ginv(i.mat.beta.s01) var.betas[i,1] <- mat.beta.s01[3*i-2,1] var.betas[i,2] <- mat.beta.s01[3*i-1,2] var.s01[i] <- mat.beta.s01[3*i,3] } # Variâncias estimadas dos parâmetros estimados pelo método numérico var_c <- matrix(0,6*nrow(sigmas5),6) var.phi0 <- rep(0,nrow(sigmas5)) var.phi1 <- rep(0,nrow(sigmas5)) var.phi2 <- rep(0,nrow(sigmas5)) var.eta <- rep(0,nrow(sigmas5)) var.nu1 <- rep(0,nrow(sigmas5)) var.nu2 <- rep(0,nrow(sigmas5)) for(i in 1:nrow(sigmas5)){ var_c[(6*i-5):(6*i),] <- ginv(hessianas5[(6*i-5):(6*i),],tol=1e-5) var.eta[i] <- var_c[6*i-5,1] var.nu1[i] <- var_c[6*i-4,2] var.nu2[i] <- var_c[6*i-3,3] var.phi0[i] <- var_c[6*i-2,4] var.phi1[i] <- var_c[6*i-1,5] var.phi2[i] <- var_c[6*i,6] } # VARIÂNCIAS ESTIMADAS DO s1, s02, s2 var.s02 <- rep(0,nrow(sigmas5)) var.s1 <- rep(0,nrow(sigmas5)) var.s2 <- rep(0,nrow(sigmas5)) for(i in 1:nrow(sigmas5)){ var.s1[i] <- (re_parametros5[i,2]^2)*var.s01[i] + (sigmas5[i,1]^2)*var.nu1[i] var.s02[i] <- (re_parametros5[i,1]^2)*var.s01[i] + (sigmas5[i,1]^2)*var.eta[i] var.s2[i] <- (re_parametros5[i,3]^2)*var.s01[i] + (sigmas5[i,1]^2)*var.nu2[i] } ############################# ANÁLISE DO BCRM ################################## ## Funções de simulação, estimação e predição com BCRM rBCRM_sim <- function(n1,n2,N,n_arm,perc_co_loc=1,medias=c(0,0),sigmas,phis, kappas=c(.5,.5),sem_grid=1,sem_cam_ale=1){ set.seed((sem_grid*3331)) grx1 <- runif(n1) set.seed((sem_grid*3771)) gry1 <- runif(n1) gr1 <- unname(cbind(grx1[1:n1],gry1[1:n1])) set.seed((sem_grid*3713)) grx2 <- runif(n2) set.seed((sem_grid*3317)) gry2 <- runif(n2) gr2 <- unname(cbind(grx2,gry2)) set.seed((sem_grid*3137)) grx3 <- runif(n_arm) set.seed((sem_grid*3337)) gry3 <- runif(n_arm) gr3 <- unname(cbind(grx3,gry3)) n <- n1+n2 n_co <- round(min(n1,n2)*perc_co_loc,0) k1 <- c((n1+n_co+1):n) k2 <- c(1:(n2-n_co)) k1_1 <- c((n_co+1):n1) k2_1 <- c(1:(n1-n_co)) k3 <- c(1:n) k4 <- c(1:n1) k4_1 <- c(1:(n1+n_arm)) k5 <- c((n1+1):n) k5_1 <- c((n1+n_arm+1):(n+2*n_arm)) k5_2 <- c((n1+n_arm+1):(n1+n_arm+n2)) k6 <- c(1:n2) k6_1 <- c(1:(n2+n_arm)) k7 <- c((n1+n_arm+n_co+1):(n+n_arm)) k8 <- c(1:(n+2*n_arm)) if(n1==n2) gr <- rbind(gr1,gr3,gr1,gr3) if(n1==n2 & perc_co_loc < 1) gr[k7,] <- gr2[k2,] if(n1>n2) gr <- rbind(gr1,gr3,gr1[k6,],gr3) if(n1>n2 & perc_co_loc < 1) gr[k7,] <- gr2[k2,] if(n2>n1) gr <- rbind(gr2[k4,],gr3,gr2,gr3) if(n2>n1 & perc_co_loc < 1) gr[k1_1,] <- gr1[k2_1,] m_1 <- medias[1] m_2 <- medias[2] s11 <- sigmas[1] s21 <- sigmas[2] s22 <- sigmas[3] phi1 <- phis[1] phi2 <- phis[2] kappa1 <- kappas[1] kappa2 <- kappas[2] set.seed(sem_cam_ale*31) R11 <- grf(n=(n1+n_arm),grid=gr[k4_1,],nsim=N,cov.pars=c(s11^2,phi1),kappa=kappa1) set.seed(sem_cam_ale*31) R21 <- grf(n=(n2+n_arm),grid=gr[k5_1,],nsim=N,cov.pars=c(s21^2,phi1),kappa=kappa1) set.seed(sem_cam_ale*37) R22 <- grf(n=(n2+n_arm),grid=gr[k5_1,],nsim=N,cov.pars=c(s22^2,phi2),kappa=kappa2) y1 <- matrix(0,(n1+n_arm),N) y2 <- matrix(0,(n2+n_arm),N) if(N >1){ y1 <- m_1 + R11$data y2 <- m_2 + R21$data + R22$data } if(N == 1){ y1[,1] <- m_1 + R11$data y2[,1] <- m_2 + R21$data + R22$data } saida <- list() geo1 <- cbind(gr[k4,],y1[k4,]) geo2 <- cbind(gr[k5_2,],y2[k6,]) geo_arm1 <- cbind(gr3,y1[(n1+1):(n1+n_arm),]) geo_arm2 <- cbind(gr3,y2[(n2+1):(n2+n_arm),]) saida$geo_arm1 <- as.geodata(geo_arm1,data.col=c(3:(2+N))) saida$geo_arm2 <- as.geodata(geo_arm2,data.col=c(3:(2+N))) saida$geo1 <- as.geodata(geo1,data.col=c(3:(2+N))) saida$geo2 <- as.geodata(geo2,data.col=c(3:(2+N))) saida$geo1$coords <- unname(saida$geo1$coords) saida$geo2$coords <- unname(saida$geo2$coords) saida$geo_arm1$coords <- unname(saida$geo_arm1$coords) saida$geo_arm2$coords <- unname(saida$geo_arm2$coords) saida$teste1 <- R11$data saida$teste2 <- R21$data return(saida) } emvBCRM_sim <- function(geo1,geo2,corr_ini,ruido=F,kappas=c(0.5,0.5),metodo="Nelder-Mead",inferior=-Inf,superior=Inf, controle=list(),hessiana=F){ n1 <- nrow(geo1$coords) n2 <- nrow(geo2$coords) n <- n1+n2 N <- ncol(geo1$data) if(is.null(N)) N <- 1 k1 <- c(1:n1) k2 <- c((n1+1):n) k3 <- c(1:n) k4 <- c(1:n2) dados <- matrix(0,n,N) if(N > 1) for(i in 1:N) dados[,i] <- c(geo1$data[,i],geo2$data[,i]) if(N == 1) dados[,1] <- c(geo1$data,geo2$data) dists <- unname(as.matrix(dist(rbind(geo1$coords,geo2$coords),diag=TRUE,upper=TRUE))) X <- matrix(c(rep(1,n1),rep(0,n2),rep(0,n1),rep(1,n2)),n,2) kappa1 <- kappas[1] kappa2 <- kappas[2] ctes <- list() ctes$n1 <- n1 ctes$n2 <- n2 ctes$dados <- dados ctes$dists <- dists ctes$X <- X est_num<-function(p,Ct,rui=F,kap,repet){ if ( any(p < 1e-20)) return(.Machine$double.xmax^.5) if ( any(p[3:4] > 1)) return(.Machine$double.xmax^.5) if ( any(p > 100)) return(.Machine$double.xmax^.5) R1 <- (((2^(kap[1]-1))*gamma(kap[1]))^(-1)) * ((Ct$dists/p[3])^(kap[1])) * besselK(Ct$dists/p[3],nu=kap[1]) R1[ Ct$dists < 0.005 ] <- 1 R2 <- (((2^(kap[2]-1))*gamma(kap[2]))^(-1)) * ((Ct$dists/p[4])^(kap[2])) * besselK(Ct$dists/p[4],nu=kap[2]) R2[ Ct$dists < 0.005 ] <- 1 V1 <- R1[k1,k1] V2 <- p[1]^2 * R1[k2,k2] + p[2]^2 * R2[k2,k2] V3 <- p[1] * R1 [k1,k2] V <- cbind(rbind(V1,t(V3)),rbind(V3,V2)) if(rui == TRUE) {diag(V) <- diag(V) + c(rep((p[7]^2),Ct$n1),rep((p[8]^2),Ct$n2))} Vinv <- ginv(V) tVinvX <- crossprod(Vinv,Ct$X) Q_hat <- drop(sum(crossprod(Ct$dados[,repet],Vinv)*Ct$dados[,repet]) - crossprod(Ct$dados[,repet],tVinvX)%*%ginv(crossprod(tVinvX,Ct$X))%*%crossprod(tVinvX,Ct$dados[,repet])) lQ_hat <- log(Q_hat) ldv <- log(det(V)) if(is.infinite(lQ_hat) | is.nan(lQ_hat) | is.infinite(ldv) | is.nan(ldv)) return(.Machine$double.xmax^.5) ll <- drop(-0.5*((Ct$n1+Ct$n2)*(log(2*pi)+lQ_hat-log(Ct$n1+Ct$n2)+1)+ldv)) if(is.infinite(ll) | is.nan(ll)) return(.Machine$double.xmax^.5) return(-ll) } resultados <- list() resultados$betas <- matrix(0,N,2) colnames(resultados$betas) <- c("mu1","mu2") resultados$sigmas <- matrix(0,N,3) colnames(resultados$sigmas) <- c("s11","s21","s22") resultados$phis <- matrix(0,N,2) colnames(resultados$phis) <- c("phi1","phi2") if(ruido == TRUE){ resultados$taus <- matrix(0,N,2) colnames(resultados$taus) <- c("tau1","tau2") resultados$re_parametros <- matrix(0,N,4) colnames(resultados$re_parametros) <- c("nu1","nu2","psi1","psi2") if(hessiana == TRUE){ resultados$hessiana <- matrix(0,(6*N),6) colnames(resultados$hessiana) <- c("nu1","nu2","phi1","phi2","psi1","psi2") rownames(resultados$hessiana) <- rep(c("nu1","nu2","phi1","phi2","psi1","psi2"),N) } } if(ruido==FALSE){ resultados$re_parametros <- matrix(0,N,2) colnames(resultados$re_parametros) <- c("nu1","nu2") if(hessiana == TRUE){ resultados$hessiana <- matrix(0,(4*N),4) colnames(resultados$hessiana) <- c("nu1","nu2","phi1","phi2") rownames(resultados$hessiana) <- rep(c("nu1","nu2","phi1","phi2"),N) } } resultados$log_lik <- matrix(0,N,1) resultados$kappas <- c(kappa1=kappas[1],kappa2=kappas[2]) resultados$erros <- matrix(0,N,5) colnames(resultados$erros) <- c("repet","nu1","nu2","phi1","phi2") convergencia <- rep(0,N) for(i in 1:N){ var1_amost <- var(geo1$data[,i]) var2_amost <- var(geo2$data[,i]) var_ini <- NULL var_ini[1] <- sqrt(var1_amost) var_ini[2] <- sqrt(0.8*var2_amost) var_ini[3] <- sqrt(0.2*var2_amost) repar_ini <- 0 repar_ini[1] <- var_ini[2]/var_ini[1] repar_ini[2] <- var_ini[3]/var_ini[1] repar_ini[3] <- corr_ini[1] repar_ini[4] <- corr_ini[2] if(ruido == TRUE) { repar_ini[5] <- var_ini[4]/var_ini[1] repar_ini[6] <- var_ini[5]/var_ini[1] } estimativas <- try(optim(p=repar_ini,fn=est_num,Ct=ctes,rui=ruido,kap=c(kappa1,kappa2),repet=i, method=metodo,lower=inferior,upper=superior,control=controle,hessian=hessiana)) if(class(estimativas) == "try-error"){ resultados$erros[i,] <- c(i,NaN,NaN,NaN,NaN) resultados$betas[i,] <- c(NaN,NaN) resultados$sigmas[i,] <- c(NaN,NaN,NaN) resultados$phis[i,] <- c(NaN,NaN) if(ruido == TRUE){ resultados$taus[i,] <- c(NaN,NaN) resultados$re_param[i,] <- c(NaN,NaN,NaN,NaN) } if(ruido==FALSE) resultados$re_parametros[i,] <- c(NaN,NaN) resultados$log_lik[i,] <- NaN if(hessiana == TRUE){ resultados$hessiana[((length(estimativas$par)*i)-(length(estimativas$par)-1)):(i*length(estimativas$par)),] <- NaN } } if(class(estimativas) != "try-error"){ if(estimativas$convergence == 0){ convergencia[i] <- 1 } print(i) R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/estimativas$par[3])^(kappa1)) * besselK(dists/estimativas$par[3],nu=kappa1) R1[dists < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists/estimativas$par[4])^(kappa2)) * besselK(dists/estimativas$par[4],nu=kappa2) R2[dists < 0.005] <- 1 V1 <- R1[k1,k1] V2 <- estimativas$par[1]^2 * R1[k2,k2] + estimativas$par[2]^2 * R2[k2,k2] V3 <- estimativas$par[1] * R1[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) VinvX <- ginv(V)%*%X mu <- drop(ginv(crossprod(VinvX,X))%*%crossprod(VinvX,dados[,i])) dif_y_mu <- dados[,i]-(X%*%mu) sigmasq <- drop((1/n) * crossprod(dif_y_mu,ginv(V)) %*% dif_y_mu) #if(sigmasq < 0) sigmasq <- drop((1/n) * crossprod(dif_y_mu,chol2inv(V)) %*% dif_y_mu) if(sigmasq < 0) resultados$erros[i,] <- c(i,estimativas$par[1],estimativas$par[2],estimativas$par[3],estimativas$par[4]) sigma11 <- sqrt(sigmasq) sigma21 <- estimativas$par[1] *sigma11 sigma22 <- estimativas$par[2] * sigma11 if(ruido == TRUE) { tau_1 <- estimativas$par[5]*sigma11 tau_2 <- estimativas$par[6]*sigma11 } resultados$betas[i,] <- mu resultados$sigmas[i,] <- c(sigma11,sigma21,sigma22) resultados$phis[i,] <- c(estimativas$par[3],estimativas$par[4]) if(ruido == TRUE){ resultados$taus[i,] <- c(tau_1,tau_2) resultados$re_param[i,] <- c(estimativas$par[1],estimativas$par[2],estimativas$par[5],estimativas$par[6]) } if(ruido==FALSE) resultados$re_parametros[i,] <- c(estimativas$par[1],estimativas$par[2]) resultados$log_lik[i,] <- -estimativas$value if(hessiana == TRUE){ resultados$hessiana[((length(estimativas$par)*i)-(length(estimativas$par)-1)):(i*length(estimativas$par)),] <- estimativas$hessian } } } resultados$erros <- subset(resultados$erros,resultados$erros[,"repet"]>0) print(c("convergência",convergencia)) class(resultados) <- "emvBGCCM" return(resultados) } predBCRM_sim <- function(geos,parametros){ geo1 <- geos$geo1 geo2 <- geos$geo2 geo_a1 <- geos$geo_arm1 geo_a2 <- geos$geo_arm2 if(nrow(parametros$erros)>0){ geo1$data <- geo1$data[,-c(parametros$erros[,1])] geo2$data <- geo2$data[,-c(parametros$erros[,1])] geo_a1$data <- geo_a1$data[,-c(parametros$erros[,1])] geo_a2$data <- geo_a2$data[,-c(parametros$erros[,1])] } grp <- geo_a1$coords coords1 <- geo1$coords coords2 <- geo2$coords n1 <- nrow(coords1) n2 <- nrow(coords2) n <- n1+n2 npr <- nrow(grp) nt <- npr + n k1 <- c(1:n1) k2 <- c(1:n2) k3 <- c(1:npr) k4 <- c((npr+1):(npr+n1)) k5 <- c((npr+1):(2*npr)) k6 <- c((npr+n1+1):nt) k7 <- c((n1+1):n) k8 <- c((npr+1):nt) k9 <- c(((2*npr)+1):((2*npr)+n)) N <- ncol(geo1$data) if(is.null(N)) N <- 1 dados <- matrix(0,n,N) if(N > 1) for(i in 1:N) dados[,i] <- c(geo1$data[,i],geo2$data[,i]) if(N == 1) dados[,1] <- c(geo1$data,geo2$data) gr_dados <- rbind(coords1,coords2) gr <- rbind(grp,gr_dados) dists <- matrix(0,nt,n) for(i in 1:nt){ dists[i,] <- spDistsN1(gr_dados,gr[i,]) } X <- matrix(c(rep(1,npr),rep(0,npr),rep(1,n1),rep(0,n2),rep(0,npr),rep(1,npr),rep(0,n1),rep(1,n2)),(2*npr+n),2) betas <- parametros$betas sigmas <- parametros$sigmas phis <- parametros$phis if(nrow(parametros$erros)>0){ betas <- matrix(c(parametros$betas[-c(parametros$erros[,1]),]),N,2) sigmas <- matrix(c(parametros$sigmas[-c(parametros$erros[,1]),]),N,3) phis <- matrix(c(parametros$phis[-c(parametros$erros[,1]),]),N,2) } if(!is.null(parametros$taus)) taus <- parametros$taus kappas <- parametros$kappas saida <- list() saida$pred1 <- matrix(0,npr,N) saida$pred2 <- matrix(0,npr,N) saida$var_pred1 <- matrix(0,npr,N) saida$var_pred2 <- matrix(0,npr,N) for(l in 1:N){ R1 <- (((2^(kappas[1]-1))*gamma(kappas[1]))^(-1)) * ((dists/phis[l,1])^(kappas[1])) * besselK(dists/phis[l,1],nu=kappas[1]) R1[dists < 0.005] <- 1 R2 <- (((2^(kappas[2]-1))*gamma(kappas[2]))^(-1)) * ((dists[k6,k7]/phis[l,2])^(kappas[2])) * besselK(dists[k6,k7]/phis[l,2],nu=kappas[2]) distaux <- dists[k6,k7] R2[distaux < 0.005] <- 1 Sig1 <- (sigmas[l,1]^2)*R1[k4,k1] Sig2 <- (sigmas[l,2]^2)*R1[k6,k7] + (sigmas[l,3]^2)*R2 Sig12 <- (sigmas[l,1]*sigmas[l,2])*R1[k4,k7] Sig <- cbind(rbind(Sig1,t(Sig12)),rbind(Sig12,Sig2)) if(!is.null(parametros$taus)){ diag(Sig) <- diag(Sig) + c(rep(taus[l,1],n1),rep(tau[l,2],n2)) } Sig_c1 <- (sigmas[l,1]^2)*R1[k3,k1] Sig_c12 <- (sigmas[l,1]*sigmas[l,2]) * R1[k3,k7] Sig_c <- cbind(Sig_c1,Sig_c12) m_p <- as.vector(X%*%betas[l,]) m_z_dadoy <- m_p[k3] - (Sig_c%*%ginv(Sig)%*%(m_p[k9]-dados[,l])) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred1[,l] <- m_z_dadoy saida$var_pred1[,l] <- (sigmas[l,1]^2) - diag(pondera_var) R2c <- (((2^(kappas[2]-1))*gamma(kappas[2]))^(-1)) * ((dists[k3,k7]/phis[l,2])^(kappas[2])) * besselK(dists[k3,k7]/phis[l,2],nu=kappas[2]) distaux <- dists[k3,k7] R2c[distaux < 0.005] <- 1 Sig_c2 <- (sigmas[l,2]^2)*R1[k3,k7] + (sigmas[l,3]^2)*R2c Sig_c21 <- (sigmas[l,1]*sigmas[l,2]) * R1[k3,k1] Sig_c <- cbind(Sig_c21,Sig_c2) m_z_dadoy <- m_p[k5] - (Sig_c%*%ginv(Sig)%*%(m_p[k9]-dados[,l])) pondera_var <- Sig_c%*%tcrossprod(ginv(Sig),Sig_c) saida$pred2[,l] <- m_z_dadoy saida$var_pred2[,l] <- ((sigmas[l,2]^2) + (sigmas[l,3]^2)) - diag(pondera_var) } return(saida) } ## Simulando sim1 <- rBCRM_sim(n1=100,n2=100,N=5,n_arm=20, perc_co_loc=1,medias=c(150,60),sigmas=c(9,5,2),phis=c(0.25,0.20), kappas=c(0.5,0.5),sem_cam_ale=3) ## Estimando est1 <- emvBCRM_sim(sim1$geo1,sim1$geo2,corr_ini=c(0.25,0.2),kappas=c(0.5,0.5), hessiana=T,metodo="Nelder-Mead",controle=list(maxit=10000)) ## Krigagem pr1 <- predBCRM_sim(geos=sim1,parametros=est1) ## Análise Final betas5 <- est1$betas sigmas5 <- est1$sigmas phis5 <- est1$phis re_parametros5 <- est1$re_parametros hessianas5 <- est1$hessianas log_liks5 <- est1$log_liks krigs1 <- pr1$pred1 krigs2 <- pr1$pred2 var_krigs1 <- pr1$var_pred1 var_krigs2 <- pr1$var_pred2 X <- matrix(c(rep(1,100),rep(0,100),rep(0,100),rep(1,100)),ncol=2) dados1 <- sim1$geo1$data dados2 <- sim1$geo2$data dados <- rbind(dados1,dados2) dists <- unname(as.matrix(dist(rbind(sim1$geo1$coords,sim1$geo2$coords),diag=TRUE,upper=TRUE))) kappa1 <- 0.5 kappa2 <- 0.5 k1 <- c(1:100) k2 <- c(101:200) var.betas <- matrix(0,nrow=nrow(betas5),ncol=2) var.s01 <- rep(0,nrow(sigmas5)) mat.beta.s01 <- matrix(0,nrow=3*nrow(betas5),3) for(i in 1:nrow(betas5)){ R1 <- (((2^(kappa1-1))*gamma(kappa1))^(-1)) * ((dists/phis5[i,1])^(kappa1)) * besselK(dists/phis5[i,1],nu=kappa1) R1[dists < 0.005] <- 1 R2 <- (((2^(kappa2-1))*gamma(kappa2))^(-1)) * ((dists/phis5[i,2])^(kappa2)) * besselK(dists/phis5[i,2],nu=kappa2) R2[dists < 0.005] <- 1 V1 <- R1[k1,k1] V2 <- re_parametros5[i,1]^2 * R1[k2,k2] + re_parametros5[i,2]^2 * R2[k2,k2] V3 <- re_parametros5[i,1] * R1[k1,k2] V <- cbind(rbind(V1,t(V3)), rbind(V3,V2)) iV <- ginv(V) Q_hat <- t(dados[,i] - X%*%betas5[i,]) %*% iV %*% (dados[,i] - X%*%betas5[i,]) i.mat.beta.s01 <- matrix(0,3,3) i.mat.beta.s01[1:2,1:2] <- (t(X)%*%iV%*%X)/sigmas5[i,1]^2 i.mat.beta.s01[3,3] <- (2*nrow(dados)^2)/Q_hat i.mat.beta.s01[1:2,3] <- (2*(t(X)%*%iV%*%dados[,i] - t(X)%*%iV%*%X%*%betas5[i,])) / sigmas5[i,1]^3 i.mat.beta.s01[3,1:2] <- t(i.mat.beta.s01[1:2,3]) mat.beta.s01[(3*i-2):(3*i),] <- ginv(i.mat.beta.s01) var.betas[i,1] <- mat.beta.s01[3*i-2,1] var.betas[i,2] <- mat.beta.s01[3*i-1,2] var.s01[i] <- mat.beta.s01[3*i,3] } mat.beta.s11<-mat.beta.s01 var.s11 <- var.s01 var_c <- matrix(0,4*nrow(sigmas5),4) var.phi1 <- rep(0,nrow(sigmas5)) var.phi2 <- rep(0,nrow(sigmas5)) var.nu1 <- rep(0,nrow(sigmas5)) var.nu2 <- rep(0,nrow(sigmas5)) for(i in 1:nrow(sigmas5)){ var_c[(4*i-3):(4*i),] <- ginv(hessianas5[(4*i-3):(4*i),],tol=1e-5) var.nu1[i] <- var_c[4*i-3,1] var.nu2[i] <- var_c[4*i-2,2] var.phi1[i] <- var_c[4*i-1,3] var.phi2[i] <- var_c[4*i,4] } var.s12 <- rep(0,nrow(sigmas5)) var.s22 <- rep(0,nrow(sigmas5)) for(i in 1:nrow(sigmas5)){ var.s12[i] <- (re_parametros5[i,1]^2)*var.s11[i] + (sigmas5[i,1]^2)*var.nu1[i] var.s22[i] <- (re_parametros5[i,2]^2)*var.s11[i] + (sigmas5[i,1]^2)*var.nu2[i] }