Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » pessoais » codigo_fonte_em_r

Essa é uma revisão anterior do documento!

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


QR Code
QR Code pessoais:codigo_fonte_em_r (generated for current page)