Why Krige returns 2d Gaussian model one len_scale?

[whyzjuhit]

Hi,
I'm trying this example with fit_variogram=True, but find the model to be fit changed to a model fitted and the parameter also changes.

import numpy as np
import gstools as gs

cond_x = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_y = [1.2, 0.6, 3.2, 4.4, 3.8]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
gridx = np.arange(0.0, 5.5, 0.1)
gridy = np.arange(0.0, 6.5, 0.1)
model2 = gs.Gaussian(
dim=2, len_scale=[1,4], anis=0.2, angles=-0.5, var=0.5, nugget=0.1
)
OK2 = gs.krige.Ordinary(model2, [cond_x, cond_y], cond_val, exact=True, fit_variogram=True)

The original model is Gaussian(dim=2, len_scale=[1,4], anis=0.2, angles=-0.5, var=0.5, nugget=0.1)
The fitted model changes to Gaussian(dim=2, var=0.262, len_scale=6.05, nugget=4.2e-06...)

In another trying, I found the len_scale is quite different.

import numpy as np
from matplotlib import pyplot as plt
import gstools as gs

cond_x = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_y = [1.2, 0.6, 3.2, 4.4, 3.8]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]
gridx = np.arange(0.0, 5.5, 0.1)
gridy = np.arange(0.0, 6.5, 0.1)
model2 = gs.Gaussian(dim=2, len_scale=[1,4], var=0.5)
OK2 = gs.krige.Ordinary(model2, [cond_x, cond_y], cond_val, exact=True, fit_variogram=True)
ff2 = OK2.structured([gridx, gridy])

fig, ax = plt.subplots()

herten_trans_skip = ff2[0]
gamma_x = gs.vario_estimate_axis(herten_trans_skip, direction="x")[0:50]
gamma_y = gs.vario_estimate_axis(herten_trans_skip, direction="y")[0:50]

x_plot = np.arange(len(gamma_x))
y_plot = np.arange(len(gamma_y))

fit_model_x = gs.Gaussian(dim=2)
fit_model_x.fit_variogram(x_plot, gamma_x)

fit_model_y = gs.Gaussian(dim=2)
fit_model_y.fit_variogram(y_plot, gamma_y)

ax.scatter(x_plot, gamma_x, label="estimated variogram in x-dir")
ax.plot(
	x_plot,
	fit_model_x.variogram(x_plot),
	linestyle="--",
	label="variogram in x-dir",
)
ax.scatter(y_plot, gamma_y, label="estimated variogram in y-dir")
ax.plot(
	y_plot,
	fit_model_y.variogram(y_plot),
	linestyle="--",
	label="variogram in y-dir",
)
ax.legend(loc="lower right")
ax.set_title("Fitting an anisotropic model")
plt.show()

the len_scale of fit_model_x is 1.06e2 and fit_model_y is 38.3, but the original is [1,4]. How can I fix it?

[MuellerSeb]

Hey there,

first, setting len_scale=[1,4], anis=0.2 is redundant, since anis will be set to 4 from the given length scales.

Second, there are to few data points to estimate a variogram. You could see that by doing the variogram estimation by hand:

import gstools as gs
from matplotlib import pyplot as plt

cond_x = [0.3, 1.9, 1.1, 3.3, 4.7]
cond_y = [1.2, 0.6, 3.2, 4.4, 3.8]
cond_val = [0.47, 0.56, 0.74, 1.47, 1.74]

bins, gamma = gs.vario_estimate((cond_x, cond_y), cond_val)

plt.scatter(bins, gamma)
plt.show()

This means, the variogram can't be fitted correctly.

Why do you expect the model to be the same when you set fit_variogram=True? This will alter the model inplace, so the parameter will most likely be different. If you add a print(model2) after the kriging initialization, you will see it changed.

Same problem occurs in the second example. Since the model is automatically fitted to the estimated variogram of the conditioning data, the anisotropy and length scales are way of, since the data doesn't carry enough information for that. Subsequently estimating these values on the kriged field results in useless parameters.

So, if you don't want, that the model is altered, don't set fit_variogram=True.

Hope that helps,
Sebastian

[whyzjuhit]

Hi, Sebastian,
Thanks for your reply.
Actually I want to use the Krige to estimate the field len_scale. So I choose a model, use Krige fit it, generate new data and estimate the len_scale along the axis. But the len_scale I get is different from what I set. According to your reply, more data should be used to Krige. I'll try it later and share the result with you as soon as I solve it.

[MuellerSeb]

I think this is not how kriging works. To use kriging you need the model to be setup already either by setting the parameters by hand or by fitting it to the estimated variogram from the conditioning data. Using the kriging field to estimate the len_scale is a chicken-and-egg problem since you would need the len_scale beforehand.

[whyzjuhit]

Oh, I understand. I should estimated variogram from the conditioning data first and then use the kriging field. I've tried by setting fit_variogram=False and fit_variogram with estimated len_scale, the data fits better. Thanks a lot.