CERENA - IST/UL

Screen-effect

Uma das fraquezas do método da indicatriz é que não há garantias de: \[P\{ z(s_0) < z \} \in [0, 1]\] Este resultado deve-se ao efeito-janela ou screen-effect, que ocorre quando o ponderador de uma amostra é fortemente penalizado por causa da presença de outra amostra muito próxima e ambas são colineares do nó a prever.

Nestes slides iremos ver um exemplo para ilustrar o screen-effect nos ponderadores de Kriging e o efeito negativo que pode produzir na probabilidade obtida.

Dados

Vamos prever a probabilidade de exceder um valor de corte (valor = 1) no nó da grid com coordenadas x = 29, y = 29 usando 13 amostras vizinhas:

zs = c(0, 0, 1,  1, 1, 0,  1, 1, 1, 0, 0, 0, 0)
cx = c(27, 30, 29, 27, 31, 30, 25, 32, 33, 30, 30, 28, 29)
cy = c(47, 24, 17, 29,  31, 30, 13,  14, 21, 12, 44, 29, 24 )
n = length(zs)

dados = data.frame(cx, cy, zs)
zs0 = data.frame(cx = 29, cy = 29)

O valor (0 ou 1) em cada amostra representa uma realização da variável indicatriz.

Mapa das amostras

Pares amostras (3,13), (4,12), (5,6) colineares com o nó a prever.

Cálculos

Vamos considerar que a dependência espacial tem forma exponêncial, sem anisotropia :

\[\gamma(h)=0 + 10* \big[1-\exp(-3*h / 15)], h>0\]

onde \(c_0 = 0\) , \(c_1 = 10\) e \(a = 15\) representam respetivamente o efeito pepita, o patamar e a amplitude do modelo.

Vamos calcular as matrizes de distância euclidiana, \(K\) e \(M\), para obter a solução do sistema de equações de Kriging : \[K^{-1}M = \lambda\]

onde os \(\lambda\) são os valores dos ponderadores das amostras, \(K\), \(M\) são as matrizes de distâncias estrutural.

Matriz de distâncias

# criar vetores coordenadas zs e zs0
xy = rbind(zs0[,1:2], dados[,1:2])

# matriz de distâncias
dij = as.matrix(dist(xy))

Matriz K

f_k = function( matd, c0, c1, a){
  
  # mat dist (zsi,zsj)
  dij_s = matd[-1,-1] 
  n = nrow(dij_s)
  
  # mod exponencial (zsi,zsj)
  gs = matrix(c0 + c1 * (1 - exp( - 3 * dij_s / a)),
              nrow = n); gs = round(gs,3)
  
  # matriz k (n+1 * n+1)
  k = rbind(gs, 1); k = cbind(k, 1)
  k [dim(k), dim(k)] = 0; k = round(k, 3)
  return(k)
}

Matriz M

f_m = function( matd, c0, c1, a){
  
  # mat dist (zs0,zs)
  dij_s0 = matd[-1,1] 
  n = nrow(matd)-1
  
  # mod exponencial (zs0,zs)
  gs0 = matrix(c0 + c1 * (1 - exp( - 3 * dij_s0 / a)),
               nrow = n); gs0 = round(gs0,3)
  
  # matriz m (n+1 * 1)
  m = c(gs0,1)
  m  = matrix(m, ncol = 1); m = round(m, 3)
  return(m)
}

Sistema de equações lineares

Calcular matrizes \(K\) e \(M\):

k = f_k (matd = dij, c0 = 0, c1 = 10, a = 15)
m = f_m (matd = dij, c0 = 0, c1 = 10, a = 15)

A solução do sistema de equações lineares de Kriging:

ik = solve(k)
w = ik %*% m

Ponderadores

O estimador impõe condição de não enviesamento:

# as.matrix(round(w[1:n],3), ncol = 1)
sum(w[1:n])
## [1] 1

mas não impõe que não haja ponderadores negativos (ou > 1):

w[w<0]
## [1] -0.002332029 -0.015113117 -0.006105471

Neste exercício há 3 ponderadores negativos.

O estimador linear de Kriging

# vetor n * 1 de lambdas
l = matrix(w[1:n], ncol = n)
# vetor n * 1 das amostras vizinhas
z = matrix(dados[,"zs"], nrow = n)
# produto escalar
v = l %*% z

A probabilidade de exceder o valor de corte é < 0:

round(v,3)
##        [,1]
## [1,] -0.019

Resultados no mapa

Id das amostras ( preto ), \(\lambda\) de valor positivo (verde) e negativo (vermelho), e probabilidade de excedência no nó (azul).

Conclusões

Uma amostra pode ter ponderador negativo quando outra amostra se encontra no seu caminho até ao nó (pontos colineares). Vimos neste exercício porquê e como pode influenciar negativamente os resultados.

Obviamente este problema não é exclusivo do Kriging normal da indicatriz, e pode ocorrer em qualquer variante do estimador de Kriging.

A solução mais simples consiste em “não mexer” nos ponderadores e substituir as probabilidades previstas por 0 (para probabilidades < 0) ou 1 (probabilidades > 1).