library(rgdal) map <- readOGR("mapas", "41mu2500gc") library(spdep) viz <- poly2nb(map) ## quantos vizinhos cada area tem nn <- card(viz) table(nn) plot(map) ce <- coordinates(map) plot(viz, ce, add=TRUE) which(nn==12) text(ce[44,1], ce[44,2], '12', col=2) ### grafo (elementos) n <- length(nn) head(nn) ii <- rep(1:n, nn) head(ii) jj <- unlist(viz[nn>0]) head(jj) ### cria grafo como matriz esparsa gg <- sparseMatrix(ii, jj, x=1) image(gg) ### vizinhança (grafo de vizinhança) vizb <- nb2listw(viz, style='B') range(eigvalb <- eigenw(vizb)) ### Seja We = vizinhança/max(autovetores) We <- gg/max(eigvalb) ### Matriz de precisao do tipo CAR C <- Diagonal(n) - 0.9*We ### Matriz de precisao do tipo SAR B <- crossprod(C) ### Algoritmo para simular de N(0, Q) ### calcule L*, tal que, L*'L* = Q (Cholesky esparso) ### simula z_1, ..., z_n, com z_i~N(0,1) ### resolva L*y = z ### Logo: Cov(y) = Cov(L*^{-1} z) = L*'Cov(z)L* = L*'L* = Q ### simula dados com cada uma (usa mesma 'semente') ### vetor N(0,1) z <- rnorm(n) ### decomposicao da matriz de precisao chol.c <- Cholesky(C) chol.b <- Cholesky(B) ### solucao do sistema linear esparso x1 <- drop(solve(chol.c, z)) x2 <- drop(solve(chol.b, z)) summary(x1) summary(x2) ### prepara para visualizar map$x1 <- x1 map$x2 <- x2 library(gridExtra) ## Observe as figuras e responda ### há correlação espacial? R.: Sim (rho usado = 0.9) ### alguma outra diferença? R.: Sim, x2 é mais suave que x1 grid.arrange(spplot(map, c('x1')), spplot(map, c('x2'))) ### teste: W <- matrix(c(0, 1/2, 1/3, 0, 1/2, 0, 1/3, 0, 1/2, 1/2, 0, 1/1, 0, 0, 1/3, 0), 4) W rowSums(W) colSums(W) tt <- c(2, 2, 3, 1) W*tt