# Ref:volta.cokri.quad.R # # carregando pacote que calcula abscissas e pesos da quadratura # de Gauss-Hermite require(statmod) volta.quad <- function(med.cov, # media e covariancia de cokrigagem n.pontos=6 # ordem da quadratura de Gauss-Hermite ){ mu.ck <- med.cov[[1]] sigma.ck <- med.cov[[2]] desvio.pd <- sqrt(diag(sigma.ck)) ic.mu <- data.frame(mu.ck-qnorm(0.975)*desvio.pd,mu.ck+qnorm(0.975)*desvio.pd) names(ic.mu) <- c("L.Minimo","L.Maximo") # Funcao que retorna uma matriz 3x4, cujo vetor de medias (resultante # de g1) sera alocado na primeira coluna(pos=1:3) e a matriz de # covariancia (resultante de g2) nas colunas de 2 a # 4(pos=4:12) g <- function(Y,mu,R,pos){ arg.agl <- mu+sqrt(2)*t(R)%*%Y g1 <- pi^-1*agl(arg.agl) g2 <- pi^-1*(agl(arg.agl)-g1)%*%t((agl(arg.agl)-g1)) resultado <- t(rbind(g1,g2)) #return(resultado[pos]) } # calculando a quadratura de gauss-hermite quad.gauss.comp <- function(mu,R,func=g,pos,np=n.pontos){ AbcPeso <- gauss.quad(np,kind='hermite') Y <- AbcPeso$nodes peso<- AbcPeso$weights soma <- 0 for(i in 1:np){ for(j in 1:np){ soma<-peso[i]*peso[j]*func(Y=c(Y[i],Y[j]),mu=mu,R=R,pos=pos)+soma } } #return(soma) } # calculando a quadratura de gauss-hermite para o vetor de medias n.linhas <- length(mu.ck)/2 compo.md <- matrix(ncol=3,nrow=n.linhas) seq1 <- seq(1,n.linhas*2,by=2) seq2 <- seq(2,n.linhas*2,by=2) for(i in 1:n.linhas){ sigma.ck <- med.cov[[2]][seq1[i]:seq2[i],seq1[i]:seq2[i]] mu.ck <- med.cov[[1]][seq1[i]:seq2[i]] cp <- quad.gauss.comp(mu=mu.ck,R=chol(sigma.ck),func=g,pos=c(1:3),np=n.pontos) compo.md[i,] <- cp } #rowSums(compo.md) # calculando a quadratura de gauss-hermite para a matriz de covariancias compo.var <- matrix(ncol=3,nrow=n.linhas) for(i in 1:n.linhas){ sigma.ck <- med.cov[[2]][seq1[i]:seq2[i],seq1[i]:seq2[i]] mu.ck <- med.cov[[1]][seq1[i]:seq2[i]] cp <- quad.gauss.comp(mu=mu.ck,R=chol(sigma.ck),func=g,pos=c(4,8,12),np=n.pontos) compo.var[i,] <- cp } retorna <- list() retorna[[1]] <- compo.md retorna[[2]] <- compo.var return(retorna) } #> simplex <- volta.quad(med.cov=md.cov.ck,n.pontos=6) #> simplex #[[1]] # [,1] [,2] [,3] # [1,] 0.2780260 0.2282011 0.4937729 # [2,] 0.2691846 0.2294145 0.5014009 # [3,] 0.2622219 0.2305321 0.5072461 # [4,] 0.2561579 0.2317903 0.5120518 # [5,] 0.2499043 0.2334437 0.5166520 # [6,] 0.2513668 0.2333344 0.5152988 # [7,] 0.2599593 0.2316183 0.5084223 # [8,] 0.2668678 0.2306520 0.5024802 # [9,] 0.2722104 0.2304115 0.4973781 # [10,] 0.2759495 0.2309633 0.4930873 # [11,] 0.2691383 0.2294203 0.5014414 # [12,] 0.2616243 0.2303668 0.5080089 # [13,] 0.2569343 0.2310930 0.5119726 # [14,] 0.2540445 0.2318098 0.5141456 # [15,] 0.2529865 0.2324410 0.5145725 # # # [97,] 0.2746740 0.2295313 0.4957947 # [98,] 0.2852018 0.2275748 0.4872234 # [99,] 0.2967117 0.2256638 0.4776246 #[100,] 0.3088924 0.2237331 0.4673745 #[[2]] # [,1] [,2] [,3] # [1,] 0.04250532 0.02620222 0.1246324 # [2,] 0.03996517 0.02651572 0.1280613 # [3,] 0.03782427 0.02679178 0.1305852 # [4,] 0.03572346 0.02708526 0.1324765 # [5,] 0.03328306 0.02745801 0.1340499 # [6,] 0.03360839 0.02742101 0.1333780 # [7,] 0.03648064 0.02701017 0.1306110 # [8,] 0.03844333 0.02675710 0.1279395 # [9,] 0.03947674 0.02664841 0.1253067 # [10,] 0.03946644 0.02669246 0.1226423 # [11,] 0.03995246 0.02651726 0.1280799 # [12,] 0.03787743 0.02676970 0.1310894 # [13,] 0.03646492 0.02695166 0.1328204 # [14,] 0.03536436 0.02711284 0.1335950 # [15,] 0.03464316 0.02723914 0.1334774 # # [98,] 0.04420545 0.02600606 0.1214431 # [99,] 0.04785008 0.02554449 0.1173562 #[100,] 0.05162834 0.02507282 0.1129311