CE-097 : 1o Semestre de 2015

Código visto em aula

require(geoR)
## conjunto da dados "bem comportado" a ser utilizado no exemplo desta aula
data(s100)
plot(s100)
plot(s100, low=T)
## outros dados com comportamentos diferentes
plot(Ksat)
plot(ca20)
plot(parana)
##
plot(s100, low=T)
points(s100)
## estimação dos parâmetros por máxima verossimilhança
ml <- likfit(s100, ini=c(var(s100$data), 0.1))
ml
## definindo uma malha de pontos para predição espacial (mapa de predição)
## uma malha "grossa"de predição
g1 <- expand.grid(seq(0,1, l=11), seq(0,1, l=11))
points(s100)
points(g1, col=2, pch=19, cex=0.3)
## uma malha "fina" de predição
g2 <- expand.grid(seq(0,1, l=51), seq(0,1, l=51))
points(s100)
points(g, col=2, pch=19, cex=0.3)
## predição na malha grossa:
kr1 <- krige.conv(s100, loc=g1, krige=krige.control(obj.m=ml))
image(kr1)
points(s100, add=T)
## predição na malha "fina"
kr2 <- krige.conv(s100, loc=g2, krige=krige.control(obj.m=ml))
image(kr2)
## pode-se usar outras palhetes de cores
image(kr2, col=gray(seq(1, 0, l=50)))
image(kr2, col=gray(seq(1, 0.2, l=50)))
points(s100, add=T)
image(kr2, col=gray(seq(0, 1, l=50)))
points(s100, add=T)
image(kr2, col=rev(gray(seq(0, 1, l=50))))
image(kr2, col=terrain.colors(21))
## mapa de erros padrão de predição
image(kr2, col=terrain.colors(21), val = sqrt(kr2$krige.var))
points(s100, add=T)

1. No exemplo da aula passada inspecione os objetos retornados pelas funções e em especial o retornado por krige.conv(). Para este último o resultado é uma lista na qual os dois primeiros elementos são e média e variãncia das distribuições condicionais (preditivas). Considere antão a obtenção dos seguintes tópicos adicionais:
   1. Obtenha um mapa de divida dos valores preditos nos intervalos [-Inf a -1,5), [-1,5 a 0], [0, a 1,5] e [1,5 a Inf].
   2. Obtenha um mapa de probabilidades de P[Y > 1], P[Y > 1.5], P[Y > 2].
   3. Obtenha um mapa do 1o quartil, um do 3o quantil e dos percentis 10 e 90.
   4. Como voce faria para obter a predição da proporção da área acima de um determinado limiar (por exemplo 1,2).
2. Explore um segundo exemplo de análises
   1. Comandos do R
   2. Arquivo de dados
3. Mais um exemplo:

   class(Ksat)
   names(Ksat)
   summary(Ksat)
   plot(Ksat)
   plot(Ksat, lam=0)
    
   par(mfrow=c(1,1))
   K.v <- variog(Ksat, max.dist=9, lam=0)
   plot(K.v)
    
   var(log(Ksat$dat))
    
   K.eye <- eyefit(K.v)
   K.eye
    
   gr <- expand.grid(seq(0, 22.5, by=0.2), seq(0,13, by=0.2))
   kr <- krige.conv(Ksat, loc=gr, krige=krige.control(cov.pars=c(2.3,2.35), lam=0))
   image(kr, col=gray(seq(1,0,l=21)))
   image(kr, col=terrain.colors(16))

Como ficariam perguntas semelhantes às feitas para o 1o exemplo neste último caso? Por exemplo: mapa da P[Y > 1] , mapa de quantis etc