Consider that we have two independent Matérn fields with different parameters θa and θband that the fields are separated from each other by a parametric curve (for example, a straight line for simplicity, but it could be more complex, that would simply just increase the amount of parameters to be estimated). These two fields could correspond to, for example, sea and land, where the curve separating them would represent a shoreline.
Let the composite field of the two fields be a two-dimensional 10 ×10 (for example, kilometers) zero mean Matérn field, which is discretised using a128×128grid resulting in16384individual points as before. Similarly as in Section 4.1, we shall fix σ = 1for both fields. Let bothν = 1.5and` = 1.5for one field, andν = 1.8and` = 0.9for the other field in the Matérn covariance function. Also let the interceptb and the slopek of the straight separating line be 12 and -1.5 respectively.
4.2.1 Drawing realisations
Realisations of this composite field can be drawn in a similar way as was described in Section 4.1.1. In Figure 12 are two different realisations of the composite random field.
The field to the left of the separating line is the field with the parameter value θa = [νa, `a] = [1.5,1.5]and the field to the right is the field with the parameter value θb = [νb, `b] = [1.8,0.9].
(a) (b)
Figure 12: Two realisations of the composite random field.
4.2.2 Parameter estimation with missing data
Similarly as in Section 4.1.2, missing measurement data can be simulated by omitting parts of the data. Figure 13 shows the realisation in Figure 12b with the same two elliptical regions of missing measurements as in Figure 4.
Because of the independence assumption, the likelihood of a realisation f of the com-posite field given parametersθa, θb, k and b, is p(f|θ) = N(fa|0,Cθa)N(fb|0,Cθb), whereθ = [k, b,θa,θb], 0are zero vectors of suitable sizes,fa is the left part andfb is the right part of the composite field realisation according to the parametersk andb, Cθa is the Matérn covariance matrix for the left part andCθb is the Matérn covariance matrix for the right part. According to Equation (3), the posterior ofθis then proportional to this likelihood multiplied by the prior of the parameterθ. Let us use a uniform distribution from 0 to 10 for ν and a uniform distribution from 0 to 10 for ` of both field parts, a uniform distribution from -3 to 1 forkand a uniform distribution from 5 to 20 forbas the prior. Given a single realisation, the parameterθ can be estimated using Gibbs sampling
with AM. Let us estimate the parameterθ using the realisation shown in Figure 13. The MAP estimate of the parameterθcan be used as the starting point for the Gibbs sampling algorithm.
Figure 13: Realisation of the composite field with missing data.
Figure 14 depicts the samples drawn from the posterior ofθ using Adaptive Metropolis-within-Gibbs in the case of missing data and a composite field. The computation of the samples took some 41 hours.
-1.505 -1.5 -1.495 Figure 14: Samples drawn from the posterior ofθwith missing data.
Figures 15, 16 and 17 describe the chain of samples drawn from the posterior ofθand the marginal histograms in the case of missing data and a composite field. The chain seems to mix just fine and the desired distribution is reached rather fast. The conditional sample
mean ofθcalculated the same way as previously, is [-1.500, 11.996, 1.495, 1.519, 1.805, 0.899].
101 2000 4000 6000 8000 10000 -1.506
101 2000 4000 6000 8000 10000 11.96
-1.51 -1.505 -1.5 -1.495 -1.49 k
11.94 11.96 11.98 12 12.02 12.04 b
Figure 15: Sampled chain and marginal histograms ofkandbin the case of missing data.
101 2000 4000 6000 8000 10000
101 2000 4000 6000 8000 10000 1.4
Figure 16: Sampled chain and marginal histograms ofθain the case of missing data.
101 2000 4000 6000 8000 10000
101 2000 4000 6000 8000 10000 1.7
Figure 17: Sampled chain and marginal histograms ofθbin the case of missing data.
4.2.3 Parameter estimation with missing and noisy data
Similarly as in Section 4.1.4, noise in the measurement data is simulated by adding i.i.d.
Gaussian noise, but this time the standard deviation of the added noise is only 0.001.
Figure 18 illustrates the realisation with missing data in Figure 13 with added i.i.d. Gaus-sian noise. The level of noise is not that high in this case, the standard deviation of the added noise is only around 0.04 percent of the maximum absolute value of the realisation, therefore the noise is not as apparently visible as it is in Figure 7.
Figure 18: Realisation of the composite field with missing data and added noise.
The noisy composite realisation fˆcan be expressed as fˆ = Af +, where f is the realisation of the underlying unobserved composite field,Ais a matrix which selects the non-obstructed points andis a vector of i.i.d. Gaussian noise of appropriate size. Prior of θ is the same as in the previous Section, but now the likelihood p(fˆ|θ) is according to Equation (1) expressed as N(fˆa|0,Ca +Cθa)N(fˆb|0,Cb +Cθb), whereCa and Cb are the covariance matrices of the added i.i.d. Gaussian noise of suitable sizes,Cθa is the Matérn covariance matrix of the underlying fieldFa with parameter valueθa and Cθb is the Matérn covariance matrix of the underlying fieldFbwith parameter valueθb. Let us estimate the parameterθ using the realisation in Figure 18. Figure 19 shows the samples drawn from the posterior of θ using Adaptive Metropolis-within-Gibbs in the case of missing and noisy data. The computation of the samples took roughly 49 hours.
-1.505 -1.5 -1.495
Figure 19: Samples drawn from the posterior ofθwith missing data and added noise.
Figures 20, 21 and 22 present the chain of samples drawn from the posterior of θ and the marginal histograms in the case of missing data and added noise. From Figure 20 it can be seen, that the marginal histogram of k looks approximately Gaussian, but the marginal histogram of b looks more like a uniform distribution rather than a Gaussian.
The conditional sample mean of θ, calculated the same way as previously, is [-1.500, 11.997, 1.501, 1.509, 1.811, 0.894].
After estimating the parameter θ either by using a MAP estimate or a conditional mean estimate, the underlying composite field realisation f can be reconstructed using the Equation (2) for the conditional distribution of f conditioned on fˆ, p(f|fˆ) = N(f|ΣATC−1fˆ,Σ), whereΣ= (Cθ−1+ATC−1AT)−1, using the estimated param-eterθvalue and assuming that the standard deviation of the added noise is known or has been estimated. Figure 23 presents two different reconstructions of the underlying field realisationf.
101 2000 4000 6000 8000 10000
101 2000 4000 6000 8000 10000 11.96
-1.51 -1.505 -1.5 -1.495 -1.49 k
11.94 11.96 11.98 12 12.02 12.04 b
Figure 20: Sampled chain and marginal histograms ofk andbin the case of missing data and added noise.
101 2000 4000 6000 8000 10000
101 2000 4000 6000 8000 10000 1.4
Figure 21: Sampled chain and marginal histograms ofθa in the case of missing data and added noise.
101 2000 4000 6000 8000 10000
101 2000 4000 6000 8000 10000 0.8
Figure 22: Sampled chain and marginal histograms ofθbin the case of missing data and added noise.
(a) (b)
Figure 23: Two reconstructions of the underlying composite field realisationf.
Figure 24a depicts the expected value of the underlying composite field realisation f. Comparing to Figure 24b, the actual underlying realisationf from Figure 12b, it can be seen that the regions with missing data due to obstructions have been reconstructed with some accuracy.
(a) (b)
Figure 24: Comparison of the expected value of the field realisation with the actual origi-nal realisation.