GP model of sensor measurements

Copyright (C) 2018 Juan Pablo Carbajal

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author: Juan Pablo Carbajal ajuanpi+dev@gmail.com Created: 2018-07-25

Contents

Dependencies

Load all dependencies of this script

pkg load signal      # for autocorrelation
pkg load statistics  # for sampling multivariate normal distribution

To install gpml in Octave 4.4.0 or later run

pkg install https://bitbucket.org/KaKiLa/gpml/downloads/gpml-4.2.0.tar.gz

pkg load gpml       # for gp models

Data and model description

The next plots show the raw data and a simple linear regression The results of the linear regression will be used to initialize the GP model in the subsequent sections.

data = dlmread ('data1.csv', ';', 1, 1);
t1   = data(:,1);
z1   = data(:,2);
N1   = length (z1);
pz1  = polyfit (t1, z1, 1);
z1_  = polyval(pz1, t1);

data = dlmread ('data2.csv', ';', 1, 1);
t2   = data(:,1);
z2   = data(:,2);
N2   = length (z2);
pz2  = polyfit (t2, z2, 1);
z2_  = polyval(pz2, t2);

Plots of the raw data and the linear fit

figure (1), clf
  hold on
  h = plot (t1, z1, '^', t2, z2, 'v');
  plot (t1, z1_, '-', 'color', get (h(1), 'color'))
  plot (t2, z2_, '-', 'color', get (h(2), 'color'))
  axis tight
  hold off
  xlabel('Time [d]')
  ylabel('Measurement [mV]')
  legend(h, {'Run 1','Run 2'})
s_gpmodel-1.png

Plots of the residuals, and their autocorrelation functions.

res1 = z1 - z1_;
res2 = z2 - z2_;
[acorr1 l1] = xcorr (res1, 'coeff');
[acorr2 l2] = xcorr (res2, 'coeff');

figure (2), clf
  subplot (2, 1, 1)
  hres = plot (t1, res1, '^;Run 1;', t2, res2, 'v;Run 2;');
  ylabel('Residual [mV]')
  axis tight

  # Autocorrelation of residuals
  subplot (2, 1, 2)
  hold on
  stem (l1, acorr1, 'color', get (hres(1), 'color'), 'filled');
  stem (l2, acorr2, 'color', get (hres(2), 'color'), 'filled');
  ylabel('Autocorrelation coeff.')
  xlabel('Lag [d]')
  axis tight
  hold off
s_gpmodel-2.png

Gaussian process (GP) model

The data generating model is:

$$ \dot{x} = \xi(t) $$ $$ \dot{y} = x(t) $$ $$ z_k = y(t_k) + \zeta_k$$ it is defined in $t \in [0,\infty)$ with initial conditions $x(0) = x_o, y(0) = y_o $

The random variables have the following distributions $$\xi(t) \sim \mathcal{N}(0,\delta(t-t^\prime)\sigma_{\xi}^2)$$ $$ \zeta_k \sim N(0,\sigma^2_{\zeta}) $$

This model has the following corresponding GP.

The state $x(t)$ is white noise, $$x(t) \sim \mathcal{N}(x_0, W_0)$$ where $W_0$ corresponds to the Wiener process covariance. This covariance corresponds to covW(i=0) in the gpml package.

The state $y(t)$ is the 1-time integrated Wiener process, $$y(t) \sim \mathcal{N}(x_0 t + y_0, W_1)$$ its covariance $W_1$ corresponds to covW(i=1) in the gpml package.

Finally, the observable $z_k := z(t=t_k)$ is modeled as $$ z(t) \sim \mathcal{N}(x_0 t + y_0, W_1 + \delta(t-t^\prime)\sigma^2_{\zeta}) $$

# initial values for hyper-parameters
hyp10 = hyp20 = struct ('cov', [], 'mean', [], 'lik', []);

Linear mean function for the observable

meanfunc = {@meanSum, {
                       @meanLinear, @meanConst
                      }
           };
hyp10.mean = pz1;
hyp20.mean = pz2;

The likelihood hyper-parameter is the intensity of the observational noise the value is $\log \left(\sigma_\zeta \right)$

likfunc = @likGauss;
hyp10.lik = hyp20.lik = 0;

Covariance function for the continuous observable $z(t)$ The value of the hyper-parameter is the logarithm of the standard deviation of the 1-time integrated Wiener process.

covfunc = {@covW, 1};
hyp10.cov = hyp20.cov = 0;
# The following is equivalent, if we fix the likelihood hyper-parameter to -Inf
#covfunc = {@covSum, {
#                     {@covW, 1}, @covNoise
#                     }
#           }
#hyp.cov = [0; 0];

args  = {@infExact, meanfunc, covfunc, likfunc};

Find maximum likelihood parameters

Note: the evaluation is done with evalc because the minimize function is too verbose.

if ~exist ('hyp1', 'var')
  printf ('Optimizing for Run 1 ...\n'); fflush (stdout);
  #hyp1  = minimize (hyp10, @gp, -1e3, args{:}, t1, z1);
  st = evalc ('hyp1 = minimize (hyp10, @gp, -1e3, args{:}, t1, z1);');
  printf ('%s', st(end-48:end)); fflush (stdout);
  nlml1 = gp (hyp1, args{:}, t1, z1);
endif

if ~exist ('hyp2', 'var')
  printf ('Optimizing for Run 2 ...\n'); fflush (stdout);
  #hyp2  = minimize (hyp20, @gp, -1e3, args{:}, t2, z2);
  st = evalc ('hyp2 = minimize (hyp20, @gp, -1e3, args{:}, t2, z2);');
  printf ('%s', st(end-48:end)); fflush (stdout);
  [nlml2, ~, post] = gp (hyp2, args{:}, t2, z2);
endif
Optimizing for Run 1 ...
Function evaluation    237;  Value 8.415618e+01
Optimizing for Run 2 ...
Function evaluation    141;  Value 7.538216e+01

Generate predictions in an extended interval

t         = linspace (0, max([t1(:); t2(:)])*1.25, 100).';
[z1_ dz2] = gp (hyp1, args{:}, t1, z1, t);
dz1_      = 1.96 * sqrt (dz2);
[z2_ dz2] = gp (hyp2, args{:}, t2, z2, t);
dz2_      = 1.96 * sqrt (dz2);

Plot the GP prediction

figure (3), clf
  hold on
  hz1 = shadowplot (t, z1_, dz1_);
  hz2 = shadowplot (t, z2_, dz2_);
  h = plot (t1, z1, '^', t2, z2, 'v');
  clr = get (h(1), 'color');
  set (hz1.line.center, 'color', clr);
  set ([hz1.line.top hz1.line.bottom], 'color', min (clr*0.8, 1),
    'linewidth', 1);
  set (hz1.patch, 'facecolor', min (clr*3, 1));
  #set (hz1.patch, 'facecolor', clr, 'facealpha', 0.3); # no alpha in html files
  clr = get (h(2), 'color');
  set (hz2.line.center, 'color', clr);
  set ([hz2.line.top hz2.line.bottom], 'color', min (clr*0.8, 1),
    'linewidth', 1);
  set (hz2.patch, 'facecolor', min (clr*3, 1));
  #set (hz2.patch, 'facecolor', clr, 'facealpha', 0.3); # no alpha in html files
  xlabel('Time [d]')
  ylabel('Measurement [mV]')
  legend(h, {'Run 1','Run 2'})
  axis tight
  hold off
s_gpmodel-3.png

The estimation of the quantities of interest

printf ("Run 1:\n")
printf ("    x(0) = %.2f mV/day\n", hyp1.mean(1))
printf ("    y(0) = %.2f mV\n", hyp1.mean(2))
printf ("    std. dev. y(t) = %.2f mV\n", exp (hyp1.cov))
printf ("    std. dev. \zeta(t) = %.2f mV\n", exp (hyp1.lik))
printf ("Run 2:\n")
printf ("    x(0) = %.2f mV/day\n", hyp2.mean(1))
printf ("    y(0) = %.2f mV\n", hyp2.mean(2))
printf ("    std. dev. y(t) = %.2f mV\n", exp (hyp2.cov))
printf ("    std. dev. \zeta(t) = %.2f mV\n", exp (hyp2.lik))
fflush (stdout);
Run 1:
    x(0) = -0.35 mV/day
    y(0) = 3.31 mV
    std. dev. y(t) = 0.01 mV
    std. dev. zeta(t) = 1.19 mV
Run 2:
    x(0) = -0.42 mV/day
    y(0) = 4.10 mV
    std. dev. y(t) = 0.01 mV
    std. dev. zeta(t) = 0.90 mV

Plots of the residuals, and their autocorrelation functions.

res1 = z1 - gp (hyp1, args{:}, t1, z1, t1);
res2 = z2 - gp (hyp2, args{:}, t2, z2, t2);
[acorr1 l1] = xcorr (res1);
[acorr2 l2] = xcorr (res2);

figure (4), clf
  subplot (2, 1, 1)
  hres = plot (t1, res1, '^;Run 1;', t2, res2, 'v;Run 2;');
  ylabel('Residual [mV]')
  axis tight

  # Autocorrelation of residuals
  subplot (2, 1, 2)
  hold on
  stem (l1, acorr1, 'color', get (hres(1), 'color'), 'filled');
  stem (l2, acorr2, 'color', get (hres(2), 'color'), 'filled');
  xlabel('Time [d]')
  axis tight
  hold off
s_gpmodel-4.png

Posterior (or conditioned) covariance matrices

Note: can be computed more efficiently from the structure post

function K = posteriorcov (covfunc, hyp, obs, tst)
  Koo     = feval (covfunc{:}, hyp.cov, obs);
  Kto     = feval (covfunc{:}, hyp.cov, tst, obs);
  Ktt     = feval (covfunc{:}, hyp.cov, tst);
  Koo_inv = cholinv (Koo + exp (2 * hyp.lik) * eye (length (obs)));
  K       = Ktt - Kto * Koo_inv * Kto.'; # Schur complement
endfunction

K1 = posteriorcov (covfunc, hyp1, t1, t);
K2 = posteriorcov (covfunc, hyp2, t2, t);

# Sample the posteriors distributions
z1_samp = mvnrnd (z1_, K1, 25).';
z2_samp = mvnrnd (z2_, K2, 25).';

Plot posterior mean, 95% CI, and samples from the posterior distribution

figure (5), clf
  hold on
  plot (t, [z1_samp z2_samp], '-', 'color', [0.8 0.8 0.8], 'linewidth', 1);
  hz1 = plot (t, z1_ + [0 1 -1] .* dz1_, '-', 'linewidth', 1);
  hz2 = plot (t, z2_ + [0 1 -1] .* dz2_, '-', 'linewidth', 1);
  h = plot (t1, z1, '^', t2, z2, 'v');
  clr = get (h(1), 'color');
  set (hz1(1), 'color', clr);
  set (hz1(2:3), 'color', clr * 0.8);
  set (hz1(2), 'color', get (h(1), 'color'));
  clr = get (h(2), 'color');
  set (hz2(1), 'color', clr);
  set (hz2(2:3), 'color', clr * 0.8);
  xlabel('Time [d]')
  ylabel('Measurement [mV]')
  legend(h, {'Run 1','Run 2'})
  axis tight
  hold off
s_gpmodel-5.png