Same grid

Here we demonstrate a simple case of two random fields uniformely correlated and estimation using the multivariate version of the Debiased Whittle likelihood

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 Matern32Model
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, Estimator
from debiased_spatial_whittle.grids.base import RectangularGrid
from debiased_spatial_whittle.sampling.simulation import (
    MultivariateSamplerOnRectangularGrid,
)

from debiased_spatial_whittle.models.univariate import Matern32Model 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, Estimator from debiased_spatial_whittle.grids.base import RectangularGrid from debiased_spatial_whittle.sampling.simulation import ( MultivariateSamplerOnRectangularGrid, )

Grid specification¶

In [5]:

# Note that we use the argument nvars=2 to specify that we observe a bivariate random field
g = RectangularGrid((128, 128), nvars=2)

# Note that we use the argument nvars=2 to specify that we observe a bivariate random field g = RectangularGrid((128, 128), nvars=2)

Model specification¶

we set the correlation to 0.8

In [6]:

m = Matern32Model(rho=8, sigma=1)
bvm = BivariateUniformCorrelation(m, r=0.8, f=1.5)
bvm

m = Matern32Model(rho=8, sigma=1) bvm = BivariateUniformCorrelation(m, r=0.8, f=1.5) bvm

Out[6]:

BivariateUniformCorrelation()

Name	Value	Type	Range
f	1.5	ModelParameter	(0.01, 100.0)
r	0.8	ModelParameter	(-0.99, 0.99)

Matern32Model()

Name	Value	Type	Range
rho	8	ModelParameter	(0, inf)
sigma	1	ModelParameter	(0, inf)

Sample generation¶

In [7]:

s = MultivariateSamplerOnRectangularGrid(bvm, g, p=2)
data = s()

s = MultivariateSamplerOnRectangularGrid(bvm, g, p=2) data = s()

In [8]:

fig = plt.figure()
ax = fig.add_subplot(1, 2, 1)
ax.imshow(data[..., 0], cmap="inferno")
ax = fig.add_subplot(1, 2, 2)
ax.imshow(data[..., 1], cmap="inferno")
plt.show()

fig = plt.figure() ax = fig.add_subplot(1, 2, 1) ax.imshow(data[..., 0], cmap="inferno") ax = fig.add_subplot(1, 2, 2) ax.imshow(data[..., 1], cmap="inferno") plt.show()

Profile likelihood plot for the correlation parameter¶

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

Inference¶

In [14]:

e = Estimator(db)
m = Matern32Model(rho=1, sigma=1)
m.set_param_bounds(dict(rho=(0.1, 100), sigma=(0.1, 10)))
bvm = BivariateUniformCorrelation(m, r=0.0, f=1.0)
e(bvm, data)

e = Estimator(db) m = Matern32Model(rho=1, sigma=1) m.set_param_bounds(dict(rho=(0.1, 100), sigma=(0.1, 10))) bvm = BivariateUniformCorrelation(m, r=0.0, f=1.0) e(bvm, data)

Out[14]:

BivariateUniformCorrelation()

Name	Value	Type	Range
f	1.521382992927081	ModelParameter	(0.01, 100.0)
r	0.8108108041123128	ModelParameter	(-0.99, 0.99)

Matern32Model()

Name	Value	Type	Range
rho	8.341655154633308	ModelParameter	(0.1, 100)
sigma	1.0685513540830631	ModelParameter	(0.1, 10)