Non-overlapping grids
In this example, the two random fields sampling locations do not overlap. Yet we are able to estimate the correlation parameter. Note that we would not be able to do so if the lengthscale parameter were close to zero.
Imports¶
In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
In [2]:
Copied!
from debiased_spatial_whittle.backend import BackendManager
from debiased_spatial_whittle.backend import BackendManager
In [3]:
Copied!
BackendManager.set_backend("numpy")
BackendManager.set_backend("numpy")
In [4]:
Copied!
from debiased_spatial_whittle.models.univariate import SquaredExponentialModel
from debiased_spatial_whittle.models.bivariate import BivariateUniformCorrelation
from debiased_spatial_whittle.inference.multivariate_periodogram import Periodogram
from debiased_spatial_whittle.inference.periodogram import ExpectedPeriodogram
from debiased_spatial_whittle.inference.likelihood import MultivariateDebiasedWhittle
from debiased_spatial_whittle.grids.base import RectangularGrid
from debiased_spatial_whittle.sampling.simulation import SamplerBUCOnRectangularGrid
from debiased_spatial_whittle.models.univariate import SquaredExponentialModel
from debiased_spatial_whittle.models.bivariate import BivariateUniformCorrelation
from debiased_spatial_whittle.inference.multivariate_periodogram import Periodogram
from debiased_spatial_whittle.inference.periodogram import ExpectedPeriodogram
from debiased_spatial_whittle.inference.likelihood import MultivariateDebiasedWhittle
from debiased_spatial_whittle.grids.base import RectangularGrid
from debiased_spatial_whittle.sampling.simulation import SamplerBUCOnRectangularGrid
Grid specification¶
In [5]:
Copied!
g = RectangularGrid((256, 256), nvars=2)
g.mask = np.random.rand(*g.mask.shape) > 0.2
x, y = g.grid_points
g.mask[..., 0] = np.mod(x, 40) <= 20
g.mask[..., 1] = np.mod(x, 40) > 20
g = RectangularGrid((256, 256), nvars=2)
g.mask = np.random.rand(*g.mask.shape) > 0.2
x, y = g.grid_points
g.mask[..., 0] = np.mod(x, 40) <= 20
g.mask[..., 1] = np.mod(x, 40) > 20
Model definition¶
In [6]:
Copied!
m = SquaredExponentialModel(rho=12.0)
bvm = BivariateUniformCorrelation(m, r=0.5, f=1.9)
bvm
m = SquaredExponentialModel(rho=12.0)
bvm = BivariateUniformCorrelation(m, r=0.5, f=1.9)
bvm
Out[6]:
BivariateUniformCorrelation()
| Name | Value | Type | Range |
|---|---|---|---|
f | 1.9 | ModelParameter | (0.01, 100.0) |
r | 0.5 | ModelParameter | (-0.99, 0.99) |
SquaredExponentialModel()
| Name | Value | Type | Range |
|---|---|---|---|
rho | 12.0 | ModelParameter | (0, inf) |
sigma | 1.0 | ModelParameter | (0, inf) |
Sample generation¶
In [7]:
Copied!
s = SamplerBUCOnRectangularGrid(bvm, g)
data = s()
s = SamplerBUCOnRectangularGrid(bvm, g)
data = s()
In [8]:
Copied!
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1)
ax.imshow(data[..., 0], cmap="Spectral")
ax = fig.add_subplot(1, 2, 2)
ax.imshow(data[..., 1], cmap="Spectral")
plt.show()
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1)
ax.imshow(data[..., 0], cmap="Spectral")
ax = fig.add_subplot(1, 2, 2)
ax.imshow(data[..., 1], cmap="Spectral")
plt.show()
Inference¶
In [9]:
Copied!
p = Periodogram()
p.fold = True
p = Periodogram()
p.fold = True
In [10]:
Copied!
ep = ExpectedPeriodogram(g, p)
db = MultivariateDebiasedWhittle(p, ep)
ep = ExpectedPeriodogram(g, p)
db = MultivariateDebiasedWhittle(p, ep)
In [11]:
Copied!
rs = np.linspace(-0.95, 0.95, 100)
lkhs = np.zeros_like(rs)
rs = np.linspace(-0.95, 0.95, 100)
lkhs = np.zeros_like(rs)
In [12]:
Copied!
for i, r in enumerate(rs):
bvm.r = r
lkhs[i] = db(data, bvm)
for i, r in enumerate(rs):
bvm.r = r
lkhs[i] = db(data, bvm)
In [13]:
Copied!
plt.figure()
plt.plot(rs, lkhs, "-")
plt.xlabel("correlation coefficient")
plt.ylabel("profile negative log-likelihood")
plt.show()
plt.figure()
plt.plot(rs, lkhs, "-")
plt.xlabel("correlation coefficient")
plt.ylabel("profile negative log-likelihood")
plt.show()