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> | ||
