Fits models to an experimental variogram (S.A. Harding, D.A. Murray & R. Webster).
Options
PRINT = string tokens |
Controls printed output from the fit (model, summary, estimates, correlations, fittedvalues, monitoring); default mode, summ, esti |
|---|---|
MODELTYPE = string token |
Defines which model to fit (power, boundedlinear, circular, spherical, doublespherical, pentaspherical, exponential, besselk1, gaussian, affinepower, linear, cubic, stable, cardinalsine, matern); default powe |
WEIGHTING = string token |
Method to be used for weighting (counts, cbyvar, equal); default coun |
CONSTANT = string token |
How to treat the constant (estimate, omit); default esti |
SMOOTHNESS = scalar |
Value of power parameter for the stable model, or ν parameter for the Matern model; default * i.e. estimate |
ISOTROPY = string token |
Defines whether to fit an isotropic or geometrical anisotropic model (isotropic, geometrical); default isot |
WINDOW = scalar |
Window in which to plot a graph; default 0 i.e. no graph |
TITLE = text |
Title for the graph |
XUPPER = scalar |
Upper limit for the x-axis in the graph |
PENDATA = scalar |
Pen to be used to plot the data; default 1 |
PENMODEL = scalar |
Pen to be used to plot the model; default 2 |
Parameters
VARIOGRAM = variates or matrices |
Experimental variogram to which the model is to be fitted, as a variate if in only one direction or as a matrix if there are several |
|---|---|
COUNTS = variates or matrices |
Counts for the points in each variogram (not required if WEIGHTING=equal) |
DISTANCE = variates or matrices |
Mean lag distances for the points in each variogram |
DIRECTION = variates |
Directions in which each variogram was computed |
INITIAL = scalars or variates |
Scalar defining initial distance parameter for an isotropic model, or variate with two values for a double-spherical isotropic model, or a variate with three values for a geometrical anisotropic model |
ESTIMATES = variates |
Estimated parameter values |
FITTEDVALUES = variates |
Fitted values |
EXIT = scalars |
Exit status from the nonlinear fitting |
SAVE = pointers |
Saves the model name and estimates in a pointer that can be used in KRIGE |
Description
Procedure MVARIOGRAM uses the directives FIT, FITCURVE and FITNONLINEAR to fit various models to the experimental variogram. Models must be authorized in the sense that they cannot give rise to negative variances when data are combined. Technically they are conditionally negative semi-definite (CNSD); see Webster & Oliver (1990, 2001), Journel & Huijbregts (1978) or Chiles & Delfiner (1999) for an explanation.
The MODELTYPE option specifies the model that is to be fitted. There are bounded isotropic models with finite ranges. These all take the value c + c0 for h ≥ a, and the following values for h < a
boundedlinear |
c0 + ch/a |
|---|---|
circular |
c0 + c {1 – (2/π)arccos(h/a) |
| + (2h/(πa))√(1-h2/a2)} | |
spherical |
c0 + c {1.5h/a – 0.5(h/a)3 } |
doublespherical |
c0 + c1 {1.5h/a1 – 0.5(h/a1)3 } |
| + c2 {1.5h/a2 – 0.5(h/a2)3 } for h ≤ a1 | |
| c0 + c1 + c2 {1.5h/a2 – 0.5(h/a2)3} for a1 < h < a2 | |
| where c = c1 + c2 | |
pentaspherical |
c0 + c {1.875h/a – 1.25(h/a)3 + 0.375(h/a)5} |
cubic |
c0 + c {7(h/a)2 – 8.75(h/a)3 + 3.5(h/a)5 |
| – 0.75(h/a)7} |
There are also bounded asymptotic models
besselk1 |
c0 + c {1 – h/a K1(h/a) } |
|---|---|
| (Whittle’s elementary correlation, Whittle 1954) | |
exponential |
c0 + c {1 – exp(-h/a)} |
gaussian |
c0 + c {1 – exp(-h2/a2)} |
stable |
c0 + c {1 – exp(-(h/a)b))} |
matern |
c0 + c {1 – 1 / (2(ν-1) Γ(ν)) (h/a)ν Kν(h/a)} |
unbounded models
power |
c0 + g hα |
|---|---|
| (power function with exponent α strictly between 0 and 2) | |
linear |
c0 + c h |
| which is a special case of the power function with exponent 1 |
and hole effect models
cardinalsine |
c0 + c × (1 – a/h × sin(h/a)) . |
|---|
Geometrically anisotropic models, i.e. ones that might be made isotropic by a simple linear transformation of the spatial coordinates, can be fitted by setting option ISOTROPY=geometrical. The following transformation is used:
Ω(θ) = √{ a2cos2(θ-φ) + b2sin2(θ-φ) }
where θ represents the direction (specified by the DIRECTION parameter) converted from degrees to radians. So, for example, a geometrical anisotropic power model would be
c0 + ( Ω(θ) h )power
(Note: this particular model can also be defined by setting MODELTYPE=affinepower; the ISOTROPY option is then ignored.)
In all these models, the intercept term (or nugget variance) c0 can be omitted by setting the CONSTANT option to omit; the default is estimate.
For the stable model the SMOOTHNESS option controls the power parameter for the model, and for the matern model it specifies the ν parameter. By default, the parameter is estimated. However, you can supply a value, to fix the parameter for the model fitting.
The data for the procedure can be taken directly from the FVARIOGRAM directive, with parameters DISTANCES, VARIOGRAMS and COUNTS corresponding to those with the same names in FVARIOGRAM. The data will be in variates if the variogram was calculated in only one direction. If it is in several, they can either be in matrices (as generated by FVARIOGRAM) or in variates. For ISOTROPY=geometrical, directions must be supplied, using the DIRECTIONS parameter. These should be in a variate with one value for each column if the other data are in matrices; alternatively, they should be in a variate of the same length as the other variates.
The WEIGHTING option controls the weights that are used when fitting the model. The default setting counts uses the values supplied by the COUNTS parameter, cbyvar uses the COUNTS divided by the values in VARIOGRAM, and equal uses equal weights (of one).
By default, MVARIOGRAM generates rough starting values for the parameters. If the solution does not converge there are two likely reasons. The model may be unsuited for the particular experimental variogram. For example, a bounded model is specified when the variogram is clearly unbounded, or vice versa. You should choose only models that have approximately the right shape. Alternatively, the starting values may be too far from a sensible solution. You should then supply initial values using the INITIAL parameter. For a double-spherical isotropic model, INITIAL must be set to a variate with two values representing the two distance parameters. For the other isotropic models it should be set to a scalar defining the initial distance parameter. Finally, for a geometrical anisotropic model, it should be set to a variate with three values, defining the initial values for φ, the maximum distance parameter and the minimum distance parameter.
Printed output is controlled by the PRINT option, and includes all the usual settings as in FIT, FITCURVE or FITNONLINEAR. You can also produce a high-resolution graph of the data and the fitted model, by setting the WINDOW option to the number of a suitable window. By default WINDOW is zero, and no graph is produced. The TITLE option can supply a title for the plot. Option XUPPER can define an upper value for the x-axis (i.e. distance), and PENDATA and PENMODEL can supply the numbers of the pens to be used to plot the experimental variogram and the fitted model respectively (by default 1 and 2). Alternatively, you can use the ESTIMATES parameter to save the parameter estimates, and plot the variogram and model later with the DVARIOGRAM procedure.
The SAVE parameter saves the parameter estimates and associated information required by the KRIGE directive. The FITTEDVALUES parameter saves the fitted values, and the EXIT parameter saves the exit “status code” from FIT, FITCURVE or FITNONLINEAR. A zero value indicates success (see the Guide to the Genstat Command Language, Part 2, Section 3.7.4).
Options: PRINT, MODELTYPE, WEIGHTING, CONSTANT, SMOOTHNESS, ISOTROPY, WINDOW, TITLE, XUPPER, PENDATA, PENMODEL.
Parameters: VARIOGRAM, COUNTS, DISTANCE, DIRECTION, INITIAL, ESTIMATES, FITTEDVALUES, EXIT, SAVE.
Method
The model is fitted using directives FIT, FITCURVE or FITNONLINEAR as appropriate.
Action with RESTRICT
If the data variates are restricted, only the units not excluded by the restriction will be used.
References
Chiles, J-P. & Delfiner, P. (1999). Geostatistics: Modeling Spatial Uncertainty. Wiley, Chichester.
Journel, A.G. & Huijbregts, C.J. (1978). Mining Geostatistics. Academic Press, London.
Webster, R. & Oliver, M.A. (1990). Statistical Methods in Soil and Land Resource Survey. Oxford University Press.
Webster, R. & Oliver, M.A. (2001). Geostatistics for Environmental Scientists. Wiley, Chichester.
Whittle, P. (1954). On stationary processes in the plane. Biometrika, 41, 434-449.
See also
Directives: FVARIOGRAM, MCOVARIOGRAM, KRIGE.
Procedure: DVARIOGRAM.
Commands for: Spatial statistics.
Example
CAPTION 'MVARIOGRAM example',\
!t('Model a variogram levels of potassium levels at',\
'Brooms Barn Experimental Station (see Webster & Oliver,',\
'1990, Statistical Methods in Soil and Land Resource Survey,',\
'Oxford University Press, pages 267-269).'); STYLE=meta,plain
VARIATE [VALUES=0,45,90,135] Angles
MATRIX [COLUMNS=13; ROWS=Angles] LogKvar,Kcounts,Midpoints
READ [SERIAL=yes] LogKvar,Kcounts,Midpoints
* * * *
0.005988 0.008046 0.006741 0.008364
0.008062 0.009615 0.008969 0.010654
0.011550 0.011688 0.010114 0.012356
0.013897 0.013190 0.011354 0.012891
0.015261 0.016131 0.014271 0.013980
0.016509 0.017591 0.016121 0.016478
0.018149 0.017546 0.016925 0.018412
0.019388 0.018019 0.017641 0.021293
0.020312 0.017763 0.016811 0.022991
0.020957 0.017235 0.017137 0.025295
0.020331 0.014874 0.017571 0.026301
0.019182 0.014330 0.017281 0.028884 :
0 0 0 0
396 374 399 376
362 1024 375 1032
971 612 1014 620
890 859 968 875
1336 1294 1490 1336
1227 939 1407 989
1106 1489 1323 1606
1340 1307 1684 1461
1157 1144 1599 1326
1239 1150 1929 1381
1060 896 1875 1087
879 1238 1808 1522 :
* * * *
1.000000 1.414214 1.000000 1.414214
2.000000 2.425227 2.000000 2.426056
3.106801 3.605560 3.106274 3.605559
4.081167 4.395480 4.080601 4.395571
5.189774 5.451577 5.187020 5.453624
6.159935 6.554852 6.157439 6.555201
7.138072 7.469445 7.136054 7.472196
8.236724 8.515127 8.231771 8.518978
9.211962 9.551049 9.207607 9.550405
10.309932 10.421715 10.303755 10.424751
11.282172 11.394851 11.277505 11.395464
12.259038 12.353009 12.254770 12.356227 :
CALCULATE Kcounts=Kcounts*(Midpoints<11.75)
" fit various possible models "
FOR Mod='LINEAR','SPHERICAL','EXPONENTIAL'
MVARIOGRAM [MODELTYPE=#Mod; PRINT=model,summary,estimates;\
WEIGHTING=counts] LogKvar; COUNTS=Kcounts; DISTANCES=Midpoints
ENDFOR