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]:

import numpy as np
import matplotlib.pyplot as plt

import numpy as np import matplotlib.pyplot as plt

In [2]:

from debiased_spatial_whittle.backend import BackendManager

from debiased_spatial_whittle.backend import BackendManager

In [3]:

BackendManager.set_backend("numpy")

BackendManager.set_backend("numpy")

In [4]:

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]:

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]:

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]:

s = SamplerBUCOnRectangularGrid(bvm, g)
data = s()

s = SamplerBUCOnRectangularGrid(bvm, g) data = s()

In [8]:

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]:

p = Periodogram()
p.fold = True

p = Periodogram() p.fold = True

In [10]:

ep = ExpectedPeriodogram(g, p)
db = MultivariateDebiasedWhittle(p, ep)

ep = ExpectedPeriodogram(g, p) db = MultivariateDebiasedWhittle(p, ep)

In [11]:

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]:

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]:

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()