How to use the gstools.Stable function in gstools

ax.set_xticklabels(("-1", "1"))
    ax.set_title("Direction sampling")

    ax = fig.add_subplot(122)
    #    x = np.linspace(0, np.max(rad))
    x = np.linspace(0, 10 / generator.model.integral_scale)
    y = generator.model.spectral_rad_pdf(x)
    ax.plot(x, y)
    sample_in = np.sum(rad <= np.max(x))
    ax.hist(rad[rad <= np.max(x)], bins=sample_in // 50, density=True)
    ax.set_xlim([0, np.max(x)])
    ax.set_title("Radius samples shown {}/{}".format(sample_in, len(rad)))
    fig.show()


model = Stable(dim=3, alpha=1.5)
srf = SRF(model)
plot_rand_meth_samples(srf.generator)

Exponential,
            Rational,
            Stable,
            Matern,
            Linear,
            Circular,
            Spherical,
            TPLGaussian,
            TPLExponential,
            TPLStable,
        ]
        self.std_cov_models = [
            Gaussian,
            Exponential,
            Rational,
            Stable,
            Matern,
            Linear,
            Circular,
            Spherical,
        ]
        self.dims = range(1, 4)
        self.lens = [[10, 5, 2]]
        self.anis = [[0.5, 0.2]]
        self.nuggets = [0, 1]
        self.vars = [1, 2]
        self.angles = [[1, 2, 3]]

        self.gamma_x = [1.0, 3.0, 5.0, 7.0, 9.0, 11.0]
        self.gamma_y = [0.2, 0.5, 0.6, 0.8, 0.8, 0.9]

print("Anisotropy ratios:", model.anis)
print("Main length scale:", model.len_scale)
print("All length scales:", model.len_scale_vec)
print("Main integral scale:", model.integral_scale)
print("All integral scales:", model.integral_scale_vec)


###############################################################################
# Percentile scale
# ----------------
#
# Another scale characterizing the covariance model, is the percentile scale.
# It is the distance, where the normalized
# variogram reaches a certain percentage of its sill.

model = gs.Stable(dim=3, var=2.0, len_scale=10)
per_scale = model.percentile_scale(0.9)
int_scale = model.integral_scale
len_scale = model.len_scale
print("90% Percentile scale:", per_scale)
print("Integral scale:", int_scale)
print("Length scale:", len_scale)

###############################################################################
# .. note::
#
#    The nugget is neglected by the percentile scale.
#
#
# Comparison
# ----------

.. math:: I = \int_0^\infty \rho(r) dr

You can access it by:
"""
import gstools as gs

model = gs.Stable(dim=3, var=2.0, len_scale=10)
print("Main integral scale:", model.integral_scale)
print("All integral scales:", model.integral_scale_vec)


###############################################################################
# You can also specify integral length scales like the ordinary length scale,
# and len_scale/anis will be recalculated:

model = gs.Stable(dim=3, var=2.0, integral_scale=[10, 4, 2])
print("Anisotropy ratios:", model.anis)
print("Main length scale:", model.len_scale)
print("All length scales:", model.len_scale_vec)
print("Main integral scale:", model.integral_scale)
print("All integral scales:", model.integral_scale_vec)


###############################################################################
# Percentile scale
# ----------------
#
# Another scale characterizing the covariance model, is the percentile scale.
# It is the distance, where the normalized
# variogram reaches a certain percentage of its sill.

model = gs.Stable(dim=3, var=2.0, len_scale=10)

x = np.random.RandomState(19970221).rand(1000) * 100.0
y = np.random.RandomState(20011012).rand(1000) * 100.0
model = gs.Exponential(dim=2, var=2, len_scale=8)
srf = gs.SRF(model, mean=0, seed=19970221)
field = srf((x, y))

###############################################################################
# Estimate the variogram of the field with 40 bins.

bins = np.arange(40)
bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)

###############################################################################
# Fit the variogram with a stable model (no nugget fitted).

fit_model = gs.Stable(dim=2)
fit_model.fit_variogram(bin_center, gamma, nugget=False)

###############################################################################
# Plot the fitting result.

ax = fit_model.plot(x_max=40)
ax.plot(bin_center, gamma)
print(fit_model)

scale of a covariance model. We provide two common scales with the covariance
model.

Integral scale
--------------

The `integral scale `_
of a covariance model is calculated by:

.. math:: I = \int_0^\infty \rho(r) dr

You can access it by:
"""
import gstools as gs

model = gs.Stable(dim=3, var=2.0, len_scale=10)
print("Main integral scale:", model.integral_scale)
print("All integral scales:", model.integral_scale_vec)


###############################################################################
# You can also specify integral length scales like the ordinary length scale,
# and len_scale/anis will be recalculated:

model = gs.Stable(dim=3, var=2.0, integral_scale=[10, 4, 2])
print("Anisotropy ratios:", model.anis)
print("Main length scale:", model.len_scale)
print("All length scales:", model.len_scale_vec)
print("Main integral scale:", model.integral_scale)
print("All integral scales:", model.integral_scale_vec)

###############################################################################
# Estimate the variogram of the field with 40 bins and plot the result.

bins = np.arange(40)
bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
plt.scatter(bin_center, gamma, label="data")
ax = plt.gca()

###############################################################################
# Define a set of models to test.

models = {
    "gaussian": gs.Gaussian,
    "exponential": gs.Exponential,
    "matern": gs.Matern,
    "stable": gs.Stable,
    "rational": gs.Rational,
    "linear": gs.Linear,
    "circular": gs.Circular,
    "spherical": gs.Spherical,
}
scores = {}

###############################################################################
# Iterate over all models, fit their variogram and calculate the r2 score.

for model in models:
    fit_model = models[model](dim=2)
    para, pcov, r2 = fit_model.fit_variogram(bin_center, gamma, return_r2=True)
    fit_model.plot(x_max=40, ax=ax)
    scores[model] = r2