[STATSGRASS] gstat probability maps

[Edzer J. Pebesma]

Dear Wolf & Pierluigi,

I believe that it is another common misunderstanding to believe that 
indicator kriging is needed to obtain (estimates of) probabilities. You 
can of course do it, try something along the lines of

v <- variogram(I(zinc>500)~1, meuse) 
m <- vgm(.15, "Sph", 874, .1) # just filling in some values here
ik <- krige(I(zinc>500)~1, meuse, meuse.grid, m)

but then you give up any information on how far an observation is below 
or above 500, as data are transformed to 0 (below 500) and 1 (above 500) 
values. IMO this is valuable information, and although you may criticize 
the Gaussian assumption (on the log-scale) underlying my previous 
suggestion, at least it leads to "true" probabilities, instead of 
interpolated 0 and 1 values that may (and will) as well lie outside the 
interval [0,1].

I wouldn't go as far as advising never to use indicator kriging, but the 
automatic "for probabilities you need indicator kriging" clearly misses 
important essentials of geostatistics.

Both issues, the need to simulate or the need for indicator transforms 
to obtain probabilities, have been pushed by certain authors that know 
better, and re-arise from time to time. A third and related issue is the 
suggestion that the (ordinary) kriging variance is not a measure for 
prediction error; see the following ai-geostats thread:
http://www.mail-archive.com/ai-geostats%40jrc.it/msg02890.html

Best wishes,
--
Edzer

Wolf Bergenheim wrote:
> Pierluigi,
>
> That would be indicator kriging:
>
> I use gstat, so I can't give you any specific R commands, but the
> principle is the same: You have to create a new dataset which is 1 if
> the threshold is exceeded, and 0 if not. Now your data is 0 or 1. Then
> do normal ordinary kriging on that new data column. The resulting map
> shows the probability of exceeding the threshold. In gstat you have to
> add  extra checks to see that your result is strictly [0.0 - 1.0]
>
> --Wolf
>
> On 27.03.2007 18:18, Pierluigi De Rosa wrote:
>   
>> Dear All,
>>
>> I need to prepare a probability map of exceeding a value...
>> I can make the simulation, for example I tried to make it with data meuse:
>>
>> v <- variogram(log(zinc)~1,loc= meuse) 
>> m <- vgm(.59, "Sph", 874, .04)
>> sim <- krige(formula = log(zinc)~1, locations = meuse, model = m, newdata = 
>> meuse.grid, nmax = 15, nsim = 5)
>> names(sim)
>> [1] "sim1" "sim2" "sim3" "sim4" "sim5"
>>
>> Now How can I create a probabilility map of exceding a value (for example 5)??
>>
>> Thanks
>> Pierluigi
>>
>>     
>
>   

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-stats/attachments/20070327/ae314308/attachment.html