7. Simulation

The map of means answers one question per location, and most questions a mine asks are not of that kind. What tonnage clears the cut-off? What revenue, at prices nobody knows yet? How likely is this stope to pay? These are functions of the whole field at once, and a nonlinear function of an average is not the average of the function. That is the oldest reason geostatistics simulates. This chapter covers how simulation works here, and how derived quantities are computed from it without leaving the data structures.

7.1 What a realization is

A realization is one coherent draw of the latent field at every requested location, carrying the spatial correlation, warped out to the variable’s units. predict(..., n_sim=...) stores them beside the prediction, and three properties are guaranteed rather than hoped for.

The noise is not in them. Realizations are of the ground, which is chapter 4’s doctrine. The assay scatter is integrated into the reported value and never drawn into the ensemble. Classical simulation practice conflates the two, and this package deliberately does not.
Batch invariance. A location’s simulated value does not depend on which batch it was computed in, so predicting twice, or in pieces, gives the same ensemble.
Reproducibility. One call to geoml.set_seed(...) before the model is built fixes everything: initialization, training draws and the simulation stream. A saved model keeps the seed it drew, so its simulations replay on reload (chapter 11).

For tight simulation budgets there is one option worth knowing. GPOptions(qmc_simulations=True) swaps the random draws for quasi-random Sobol ones. Measured on Walker Lake, the ensemble mean lands 7–37× closer to the exact posterior and the tail quantiles about 25% closer at the same n_sim, and the ensemble’s own spread pays nothing for it. It is off by default only because saved models must replay the streams they recorded.

Storage respects scale. Realizations live in a chunked store that is read in bands of rows, and variable.simulation(i) reads a single realization as a column. A block model’s ensemble can run to hundreds of gigabytes, and nothing in the package, or in this manual, materializes one whole.

7.2 Derived variables: arithmetic on the ensemble

A net smelter return, a metal content, a revenue: these are functions of the simulated variables, applied realization by realization, because the answer’s uncertainty comes from the inputs’ and the function bends. The container does this itself.

The model is chapter 4’s, unchanged.

import geoml
import numpy as np

geoml.set_seed(1234)
walker, walker_grid = geoml.datasets.walker()

warping = geoml.warping.ChainedWarping(
    geoml.warping.Softplus(1),
    geoml.warping.ZScore(1),
    geoml.warping.Spline(1, knots_per_arm=4))

experts = geoml.data.inducing.grid_experts(walker_grid, 10.0, block=8)

root = geoml.latent.BasicInput(
    experts,
    transform=geoml.transform.Isotropic(50))

gp = geoml.latent.BasicGP(
    root,
    size=1,
    kernel=geoml.kernels.Spherical())

model = geoml.models.VGPNetwork(
    walker, "V",
    geoml.likelihood.Gaussian(warping),
    gp,
    options=geoml.models.GPOptions(verbose=False))

model.train_full(max_iter=250)

subset = geoml.data.PointData.from_array(
    np.asarray(walker_grid.coordinates)[::20])
model.predict(subset, n_sim=50, include_noise=True)

derive walks the parents’ realizations and applies the function once per realization. The example is a payable value: only the part of the grade above a cut-off earns anything, at a price nobody has set yet. The optional simulation= keyword receives the realization’s index, which is how that external uncertainty rides along, one draw per realization.

price = np.random.default_rng(7).normal(1.0, 0.2, size=50)

payable = subset.derive(
    "payable",
    lambda v, simulation: np.maximum(v - 300.0, 0.0) * price[simulation],
    ["V"])

print("mean of derived sims:  ",
      round(float(np.mean(subset.values("payable/prediction"))), 1))
print("function of the means: ",
      round(float(np.mean(
          np.maximum(subset.values("V/prediction") - 300.0, 0.0)
          * price.mean())), 1))

The two numbers differ, and the first is the right one. Wherever the grade straddles the cut-off, some realizations pay and some do not, and applying the hinge to the mean throws that half-paying ground away. The prediction of a derived variable is the mean of its derived realizations, never the function of the parents’ predictions. The nonlinear case is the reason the class exists.

What comes back is a DerivedVariable, and it is worth knowing how it differs from the ContinuousVariable it is built on. Everything a continuous variable can do, it can do: quantiles, cut-offs, contours, grade–tonnage curves and Zarr persistence all work unchanged. What it cannot do is go back into a model. Every model refuses it as an input, because its uncertainty is inherited from its parents rather than modelled, and the package will not let that distinction blur. Two smaller differences follow from the same idea. Metadata may join the arguments as per-location constants (["V", "_metadata/density"]), and unsimulated variables are refused outright, since there is no ensemble to walk. The function is not persisted either. A reloaded derived variable is data, and rerunning the deriving script is how it is refreshed.

In the code. container.derive(names, function, arguments) in data/containers.py, and DerivedVariable in data/variables.py. The banded walk and the column reads follow the same storage discipline as everything else in the package (chapter 10).

References

Chilès, J.-P., & Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty (2nd ed.). Wiley.

7. Simulation

7.1 What a realization is

7.2 Derived variables: arithmetic on the ensemble

Further reading

References