Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
— |
pessoais:exercicio_de_simulacao_geoestatistica_--_unioeste_13_11_07 [2008/07/15 17:34] (atual) edson criada |
||
---|---|---|---|
Linha 1: | Linha 1: | ||
+ | <code> | ||
+ | ## ###################################### | ||
+ | ## | ||
+ | ## SIMULACAO DE PROCESSO GEOESTATISTICO | ||
+ | ## (Simulacao nao condicional) | ||
+ | ## | ||
+ | ## Edson Antonio Alves da Silva | ||
+ | ## (novembro, 2007) | ||
+ | ## ###################################### | ||
+ | ## Processo Gaussiano | ||
+ | ## | ||
+ | ## Modelo: Y ~ N(mu ; Sigmasq R(phi)+Tausq I) | ||
+ | ## | ||
+ | ## ###################################### | ||
+ | ## | ||
+ | ## Metodo manual usando os comando do R | ||
+ | ## | ||
+ | ## definindo os pontos amostrais | ||
+ | grid <- expand.grid(X1=1:3,X2=1:5) | ||
+ | plot(grid) | ||
+ | plot(grid, asp=1) | ||
+ | class(grid) | ||
+ | ## | ||
+ | ## definindo a matriz de distâncias (euclidianas) | ||
+ | mat.d <- round(as.matrix(dist(grid)),2) | ||
+ | mat.d | ||
+ | ## | ||
+ | ## definindo a matriz de covariancias | ||
+ | ## | ||
+ | ## ######################################## | ||
+ | ## Fcao Correlacao: exponencial (-dist/phi) | ||
+ | phi <- 3 | ||
+ | sigma.sq <- 5 | ||
+ | mat.cov <- sigma.sq*(exp(-mat.d/phi)) | ||
+ | round(mat.cov,2) | ||
+ | ## | ||
+ | ## Sorteando uma amostra normal multivariada com media mu | ||
+ | ## tomando k amostras univariadas e independentes N(0;1) | ||
+ | ## multiplicando pela matriz de covariancias e adicionando | ||
+ | ## a media (constante) | ||
+ | mu <- 100 | ||
+ | set.seed(1956) | ||
+ | Y <- rnorm(15)%*%mat.cov+mu | ||
+ | range(Y) | ||
+ | dados <- as.data.frame(cbind(grid,Y=as.numeric(round(Y,2)))) | ||
+ | dados | ||
+ | class(dados) | ||
+ | ## | ||
+ | ## convertendo os dados em um objeto geodata para uso da geoR | ||
+ | ## | ||
+ | require(geoR) | ||
+ | dados <- as.geodata(dados) | ||
+ | class(dados) | ||
+ | dados | ||
+ | ## ######################################## | ||
+ | ## Fcao Correlacao: matern | ||
+ | ## | ||
+ | kappa <- 0.5 | ||
+ | u.phi <- mat.d/phi | ||
+ | mat.cov2 <- ifelse(mat.d > 0, | ||
+ | (((2^(-(kappa-1)))/gamma(kappa))*(u.phi^kappa)*besselK(x=u.phi,nu=kappa)), | ||
+ | 1) | ||
+ | set.seed(1956) | ||
+ | Y2 <- rnorm(15)%*%mat.cov2+mu | ||
+ | dados2 <- as.data.frame(cbind(grid,Y=as.numeric(round(Y2,2)))) | ||
+ | dados2 | ||
+ | ## | ||
+ | ## ########################################################### | ||
+ | ## | ||
+ | ## Usando a funcao GRF da GeoR | ||
+ | ## | ||
+ | require(geoR) | ||
+ | args(grf) | ||
+ | ## | ||
+ | ## Simulacao SEM tendencia (media constante) em grid REGULAR | ||
+ | ## | ||
+ | ## Definindo o grid (incrementado) | ||
+ | ## | ||
+ | coord <- expand.grid(x1=(0:9)*10+5,x2=(0:9)*10+5) | ||
+ | plot(coord$x1,coord$x2) | ||
+ | A <- matrix(c(coord[,1],coord[,2]), ncol=2, byrow=FALSE) | ||
+ | pts1 <- expand.grid(seq(16,24,2),seq(36,44,2)) | ||
+ | pts2 <- expand.grid(seq(36,44,2),seq(76,84,2)) | ||
+ | pts3 <- expand.grid(seq(66,74,2),seq(6,14,2)) | ||
+ | pts4 <- expand.grid(seq(76,84,2),seq(56,64,2)) | ||
+ | B1 <- as.matrix(pts1) | ||
+ | B2 <- as.matrix(pts2) | ||
+ | B3 <- as.matrix(pts3) | ||
+ | B4 <- as.matrix(pts4) | ||
+ | area <- rbind(A,B1,B2,B3,B4) | ||
+ | plot(area[,1],area[,2],pch=20,cex=0.5) | ||
+ | ## | ||
+ | ## Parametros arbitrarios | ||
+ | ## | ||
+ | ## n=200 (dimensao do grid) | ||
+ | ## Fcao de correlacao esferica | ||
+ | ## Contribuicao = 15 | ||
+ | ## Parametro de alcance phi = 60 | ||
+ | ## Variacao de pequena escala tausq = 0 | ||
+ | ## Parametro de trasformacao lambda = 1 (sem transformacao) | ||
+ | ## Distorcao de anisotropia: | ||
+ | ## psi_A (angulo) = 0 | ||
+ | ## psi_B (razao entre eixos) = 1 | ||
+ | ## | ||
+ | dim(area) | ||
+ | set.seed(1956) | ||
+ | n=200 | ||
+ | sigma2 <- 15 | ||
+ | tausq <- 0 | ||
+ | phi <- 60 | ||
+ | Y3 <- grf(n,grid=area,cov.pars=c(sigma2,phi), nugget=tausq,cov.model='sph') | ||
+ | points(Y3,col='gray') | ||
+ | ## | ||
+ | ## | ||
+ | ## Simulacao SEM tendencia (media constante) em grid IRREGULAR | ||
+ | ## | ||
+ | set.seed(1956) | ||
+ | Y3.a <- grf(n,grid='irreg',cov.pars=c(sigma2,phi), nugget=tausq,cov.model='sph') | ||
+ | plot(Y3.a$data,Y3.a$coords[,1]) | ||
+ | points(Y3.a,col='gray') | ||
+ | ## | ||
+ | ## Simulacao COM tendencia em grid IRREGULAR | ||
+ | ## | ||
+ | ## Modelo: mu(x)= mu + beta1*coord(x1)+beta2*coord(x2) | ||
+ | ## | ||
+ | ## a) inicialmente geramos um modelo SEM tendencia | ||
+ | set.seed(1956) | ||
+ | Y3.b <- grf(n,grid='irreg',cov.pars=c(sigma2,phi), nugget=tausq,cov.model='sph') | ||
+ | names(Y3.b) | ||
+ | beta1 <- -5 | ||
+ | beta2 <- 5 | ||
+ | mu <- 10 | ||
+ | head(Y3.b$data) | ||
+ | hist(Y3.b$data) | ||
+ | set.seed(1956) | ||
+ | Y3.b$data <- (Y3.b$data +mu)+beta1*Y3.b$coords[,1]+beta2*Y3.b$coords[,2] | ||
+ | points(Y3.b) | ||
+ | plot(Y3.b$data,Y3.b$coords[,2]) | ||
+ | ## | ||
+ | ## ############################################################## | ||
+ | ## | ||
+ | ## Processos nao gaussianos | ||
+ | ## | ||
+ | ## ############################################################## | ||
+ | ## | ||
+ | ## Distribuicao log-normal | ||
+ | ## | ||
+ | Y4 <- grf(225,grid='reg', cov.pars=c(1,0.24), lambda=0) | ||
+ | hist(Y4$data) | ||
+ | ## | ||
+ | ## Distribuicao de Poisson | ||
+ | ## | ||
+ | Y5 <- grf(225, grid="reg", cov.pars=c(1, 0.25)) | ||
+ | set.seed(1956) | ||
+ | Y5$data <- rpois(225, exp(0.5 + Y5$data)) | ||
+ | head(Y5$data) | ||
+ | image(Y5, col=gray(seq(1,0.2,l=11))) | ||
+ | text(Y5$coords[,1], Y5$coords[,2], Y5$data) | ||
+ | barplot(table(Y5$data)) | ||
+ | ## | ||
+ | ## Distribuicao t-Student com 2 G.L. | ||
+ | ## | ||
+ | ## Simulacao SEM tendencia (media constante) em grid REGULAR | ||
+ | ## (nao foi utilizado a GRF da geoR) | ||
+ | ## Fcao Correlacao: exponencial (-dist/phi) | ||
+ | coord.t <- expand.grid(x1=(0:9)*10+5,x2=(0:9)*10+5) | ||
+ | dim(coord.t) | ||
+ | mat.t <- round(as.matrix(dist(coord.t)),2) | ||
+ | phi <- 3 | ||
+ | sigma.sq <- 5 | ||
+ | mat.cov.t <- sigma.sq*(exp(-mat.t/phi)) | ||
+ | mu <- 100 | ||
+ | set.seed(1956) | ||
+ | GL <- 1 | ||
+ | ## O grau de liberdade (GL) define o decaimento da curva | ||
+ | ## Valores altos de GL equivalem a norma | ||
+ | Y.t <- rt(100,GL)%*%mat.cov.t+mu | ||
+ | range(Y.t) | ||
+ | hist(Y.t) | ||
+ | dados.t <- as.data.frame(cbind(coord.t,Y=as.numeric(Y.t,2))) | ||
+ | dados.t | ||
+ | class(dados) | ||
+ | </code> |