topic How are kriging quantile maps calculated? in ArcGIS GeoStatistical Analyst Questions

How are kriging quantile maps calculated?

WilliamHuber — Fri, 27 Dec 2013 17:43:53 GMT

The help page for types of kriging surfaces seems to say that the quantile map for any quantile q (between 0 and 100%) is constructed as the sum of the prediction map and z times the standard error (SE) map where z is the qth quantile for the standard normal distribution. This in fact is not the case for ordinary kriging when the variogram is not a pure nugget, as you can check by constructing these three maps and comparing any GA quantile map (for q differing from 50%) to the prediction and SE maps. The discrepancies--which are both positive and negative and about correct on average--vary with location, prediction, and standard error, so this does not seem to be the result of an approximate calculation. What is the software really computing? What formula does it use?

(The help page echoes material from the book Using ArcGIS Geostatistical Analyst. See especially pages 262-264. I tested using no transformations of the data and specified variograms with no measurement error.)

Re: How are kriging quantile maps calculated?

EricKrause — Mon, 30 Dec 2013 17:33:25 GMT

The formula for quantiles contains a bias correction.

I'm attaching the response from a developer as a pdf. It contains the more general quantile formula for universal kriging with a transformation.

Re: How are kriging quantile maps calculated?

Anonymous User — Tue, 31 Dec 2013 21:15:16 GMT

Original User: whuber

Thank you very much, Eric.

I am trying to connect that document to what the software produces, but am failing to do so. My understanding of the "Prediction Standard Error Map" is that it contains the square root of the estimated residual variance: that is, it's the entire square root term in the "Formula for quantile map." (Cressie's notation for this term is \sigma_k(s_0).) When no transformation occurs (as in my test cases), the formula thereby looks like it should reduce to what is in Cressie: namely, the quantile map (\widehat{Z_q}) looks as if it's supposed to be the prediction map (\hat(Y)) plus a constant times the prediction SE map. *But it is not.*

The only difference I can see between what is in the document and what I have been expecting concerns the sign of x^T\mu within the square root: on one line of your document it appears with a negative sign but in the formula for the quantile map it appears with a positive sign. Isn't that negative sign a typographical error? Even if not, the discrepancies I am seeing cannot be explained by an additional term of that nature.

Allow me to explain how I know there are problems. Your formula implies W = (t(\widehat{Z_q}) - \hat{Y})^2 must be a multiple of the terms found under the square root sign. We can find that multiple by subtracting the prediction map \hat{Y} from the (transformed) quantile map t(\widehat{Z_q}) and squaring the result to produce W, and then regressing that against the squared prediction error map. I would hope to achieve a fit that is accurate at least to single precision floating point error. I haven't managed to do that--in some cases it's close but in others they are wildly different. I haven't ruled out errors on my part, but failing to find any I'm looking for all the information I can obtain concerning what the software is actually doing.

Please note that I am unconcerned with bias, accuracy, or uncertainty here: I am only trying to understand how the three GA outputs--prediction map, prediction SE map, and quantile map--are related mathematically.

Re: How are kriging quantile maps calculated?

EricKrause — Fri, 03 Jan 2014 15:16:51 GMT

Here is the response from the same developer, along with a zipped attachment of some Mathematica code, test data, and graphics:

The change in sign is not a mistyping. This is to do bias correction.

I wrote a simple test in Mathematica and the result for quantile map matches.

Attached the test:

Code.txt
Source code to simulate data with Exponential covariance function (Exponent(-3*h)).
Nugget=0
Sill=1
Range=1
Quantile=0.9
The prediction done in point (0,0) .
Remark. This code is only for testing.

OrdinaryKrigingQuantilePrediction.png �?? Wizard first page.
SemivariogramParameters.png �?? Wizard page of semivariogram parameters.
QuantilePrediction.png �?? All data are used as neighbors. Quantile prediction at point (0,0).

TestData.* - Source data in shapefile and txt formats.

Result in test code 2.15823
Result in Wizard 2.1582296949642736

If you have any more questions, it would be easier if you email me directly at ekrause@esri.com and I can put you in direct contact with the developer.