✎ pessoais:codigo_fonte_em_r

[ Por favor, preencha este local (:wiki:site_notice) ]

Você está aqui: start » pessoais » codigo_fonte_em_r

Você carregou uma revisão antiga desse documento! Se você salvá-la, irá criar uma nova versão com esses dados.

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

Por favor digite todas as letras dentro da caixa para provar que você é humano. S O K V A Por favor deixe esse campo em branco:

Resumo da edição

Observe: editando esta página você aceita disponibilizar o seu conteúdo sob a seguinte licença: GNU Free Documentation License 1.3