2d random field with mean=0 at each row of numbers and precision

[pelyme]

Hi,
i want to make a 2d random field with mean=0 at each row of numbers. But when i use gs.SRF(model, mean=0), the mean of all numbers in the random field will become zero not at each row. is there any solution for that in gstools?
And also, is there anyway to increase the precision of random field that the mean becomes 0 with less penalty?
Regards,

[MuellerSeb]

Hey there,

this is a tough one. Setting mean=0.0 actually means, that the mean at each point for an ensemble of random fields converges to this given value. Under the assumption of ergodicity, the spatial mean should be near the given mean for a large enough realization.
You are further constraining this by requiring each row to have the given mean. With increasing length of the field, this should converge.

Only thing that we are providing is a transformation to force the generated field to have the given moments (mean, variance). But this will not help with each individual row:

import gstools as gs

x = y = range(100)

model = gs.Exponential(dim=2, len_scale=5)
srf = gs.SRF(model=model, mean=0.0, seed=1234)
srf.structured([x, y])

print("mean:", srf.field.mean())
print("mean row 0:", srf.field[0].mean())

srf.transform("force_moments")

print("forced mean:", srf.field.mean())
print("forced mean row 0:", srf.field[0].mean())

Will give:

mean: 0.14957193945128852
mean row 0: -0.32238566532970436
forced mean: 1.1368683772161604e-17
forced mean row 0: -0.533112638564582

So this got even worse for the first row.

I don't see a possibility at the moment to force the given mean on each row/col of a field.

Sebastian

[pelyme]

Hi Sebastian,

Thanks for your reply and thanks for thinking about that.
By the way, is there any way for increasing the precision to get the total mean closer to the requested one?
for example in your example the mean = 0.149571939, is that possible to have a number more close to 0? (less penalty)?

Also is it possible to use different distribution such as beta distribution in generating random fields? or can I convert the data to beta distributed data?

Mohammad

[LSchueler]

Hi Mohammad,

the interesting line in Sebastians code is this one:
srf.transform("force_moments")
Try it with your own code and check the mean. Sebastian got a mean of ~1.14e-17. The exponent is important. I think that's pretty close to 0.

Regarding your second question, we provide quite a few ready to use transformations, but not for the beta distribution. If you want to add that transformation to GSTools, we would be happy about your contribution. If not, you can apply your own functions as a transformation, as documented here.

I hope that helps!
Lennart

[pelyme]

Thanks Lennart for reminding the "force_moments" part. i have missed that. really helpful.

Great. will look at the transformations options. If my problem didn't get fixed by those, then will get back to you for implementing the beta distribution.

Mohammad

[MuellerSeb]

Regarding transformation to different distributions:
We already provide some other distributions like uniform, arcsin, uquad or lognormal.

Using the uniform distribution as basis, it is really easy to convert the field to any distribution provided by scipy.stats using the ppf method (see scipy.stats.beta.ppf):

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
import gstools as gs

x = y = range(100)

model = gs.Exponential(dim=2, len_scale=5)
srf = gs.SRF(model=model, mean=0.0, seed=1234)
srf.structured([x, y])

# transform to uniform distribution on [0, 1]
srf.transform("uniform")
# transform to beta distribution with ppf from uniform distribution
# 'a' and 'b' are kwargs forwarded to the function in use
srf.transform("function", function=beta.ppf, a=2.31, b=0.627)
srf.plot()

# see histogram
p = np.linspace(0, 1)
plt.figure(2)
plt.hist(srf.field.ravel(), density=True)
plt.plot(p, beta.pdf(p, a=2.31, b=0.627))
plt.show()

[MuellerSeb]

This is worth a tutorial 😉

[pelyme]

Thank you for the answer. 👍 that is very helpful.