s_delayed_decomp

Contents

# Copyright (C) 2017 - Juan Pablo Carbajal
#
# This progrm 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 .
#
# Author: Juan Pablo Carbajal 
#

Analysis of delayed superpositions

This script demonstrates the use of the mean and covariance functions to infer the number of delays in a signal generated by delayed superposition of a given waveform,

$$ y(t) = \sum_{i=1}^N \varphi(t-\tau_i) $$

The delays are smpled from a known distribution.

In this section we define the waveform used to generate the superposition and the time vector.

s   = 0.025;
phi = @(t) exp (-t.^2 / 2 / s.^2);

nT = 200;                  # Time samples
T  = 1;                    # Duration
t  = linspace (0,T, nT).'; # Time vector

We now create the supoerposition using $N$ different delays. We do not actually impose that the delays are different, but the probability of drawing two equal delas is very low.

randg ('state', double('abraccadabra'));

nd  = 30;                       # Number of delays
d_p = {2 4};                    # Parameters for the beta distribution
d   = betarnd (d_p{:}, 1, nd);  # Delays sampled from beta dist.
pdf = @(x) betapdf (x, d_p{:}); # Beta dist. as function

y_f =@(t) sum (phi(t-d), 2);    # Delayed superposition

y = y_f(t);
figure(1)
clf
hold on
plot (t, phi(t-d), '-', 'linewidth', 1, 'color', 0.8*[1,1,1]);
line ([d;d], [zeros(1,nd); ones(1,nd)], 'color', [0.5,0.5,1])
plot (t, y,'linewidth', 2);
hold off
axis tight
xlabel ('Time')
ylabel ('y(t)')

Mean function over delays

This section implement the mean function, i.e. the average superposition over many delay arrays of the same length (the length of the array N is left as a free parameter)

$$ \bar{f}_N(t) = N \int_0^\infty \varphi(t-s) p(s) ds = N\hat{\varphi(t)} $$

The function $\hat{\varphi}$ can be evaluated on an array of observarion points, and is parametrized with the waveform and the distribution of delays.

hatphi = @(t,f,p) arrayfun (...
                 @(y)quadgk (@(s) f (y-s) .* p (s), 0, T), ...
                 t);
figure (2)
clf
h = plot (t, [phi(t-0.5) pdf(t)/max(pdf(t))], 'linewidth', 1);
set (h(1), 'color', 0.8*[1,1,1]);
set (h(2), 'color', [0.5,0.5,1]);
hold on
h(end+1) = plot (t, hatphi(t,phi,pdf),'linewidth',2);
axis tight
legend (h, {'waveform','delays pdf', 'convolution'});
xlabel ('Time')
ylabel ('a.u.')
s_delayed_decomp-1.png

Inference

In this section we implement two ways of inferring the numer of delays base don a finite set of observations and infomration about the waveform and the distribution of delays.

First, we calculate the number of delays by fittign the mean function to the observed data. We then use a Gaussian process to do the inference.

In the next section we pertrub the information to test the sensitivity of these inferences.

nTo = 15;                            # Number of observations
io  = round (linspace (1,nT,nTo)).'; # Indexes

to    = t(io);                         # Observation times
yo    = y(io);                         # Observations
sy    = 0.0;                           # Observation noise
yo    = yo + sy * randn (nTo,1);       # Noisy obs.

Provided innformation. First we consider the case of perfect information

eval (sprintf ("phi_s = %s;", func2str (phi)));
eval (sprintf ("pdf_s = %s;", func2str (pdf)));
eval (sprintf ("hatphi = %s;", func2str (hatphi)));

hphio = hatphi (to, phi_s,pdf_s);
hphi  = hatphi (t, phi_s,pdf_s);

Least squares estimation of number of delays, obtained from

$$ N_{LS} = \operatorname{argmin} \sum_{i=1}^{n} \left( y_i - N\hat{\varphi}(t)\right)^2 $$

n_ls = yo.' * hphio ./ sumsq (hphio);
y_ls = n_ls * hphi;

Gaussian process estimation of number of delays, i.e. maximum marginal likelihood solution

function [n gp] = gp_estim (f,p,t,to,io,yo,hp,hpo)
  [tt rr]   = meshgrid (to, t);
  [nT, nTo] = size (tt);

  hP = arrayfun ( ...
       @(y,x)quadgk (@(s) f (y-s) .* f (x-s) .* p (s), t(1), t(end)), ...
       tt(:), rr(:));
  hP     = reshape (hP, nT, nTo);
  hatK   = hP - hp * hpo.';
  tmp    = hatK(io,:) \ [yo hpo];
  alpha  = tmp(:,1);
  halpha = tmp(:,2);
  n      = roots ([halpha.'*hpo nTo -yo.'*alpha]);
  n      = n(n>0);

  gp     =  hatK * (alpha - n * halpha) + n * hp;
endfunction

[n_gp y_gp] = gp_estim (phi_s, pdf_s, t, to, io, yo, hphi, hphio);

function plot_result (fig, t, to, y, yo,y_ls,y_gp,nd,n_ls,n_gp)
  res = y - [y_ls y_gp];
  err = mean (res.^2);
  printf ("** MSE Error **\nLS:%.2g\t%d/%d\nGP:%.2g\t%d/%d\n", ...
          err(1), round (n_ls), nd, ...
          err(2), round (n_gp), nd);
  fflush (stdout);
  figure (fig)
  clf
  h = plot (t, [y y_ls y_gp], 'linewidth', 2);

  hold on
  h(end+1) = plot (to, yo, 'og', 'markerfacecolor', 'g');
  set (gca, 'color', 'none');
  hold off

  axis tight
  hl = legend (h, {'True', ...
              sprintf("LS Fit (n=%.1f)", n_ls), ...
              sprintf("GP Fit (n=%.1f)", n_gp), ...
              'Observed', ...
              }
         );
  set(hl, 'color', 'none');

endfunction

plot_result (3, t, to, y, yo,y_ls,y_gp,nd,n_ls,n_gp)
** MSE Error **
LS:1.4	32/30
GP:0.33	31/30
s_delayed_decomp-2.png

Repeat the results for wrongly specified pdf of delays

eval (sprintf ("phi_s = %s;", func2str (phi)));
pdf_s = @(t) normpdf (t, mean(d), std(d));
eval (sprintf ("hatphi = %s;", func2str (hatphi)));

hphio = hatphi (to, phi_s,pdf_s);
hphi  = hatphi (t, phi_s,pdf_s);

n_ls = yo.' * hphio ./ sumsq (hphio);
y_ls = n_ls * hphi;
[n_gp y_gp] = gp_estim (phi_s, pdf_s, t, to, io, yo, hphi, hphio);

plot_result (4, t, to, y, yo,y_ls,y_gp,nd,n_ls,n_gp)
** MSE Error **
LS:1.6	31/30
GP:0.33	32/30
s_delayed_decomp-3.png

Repeat the results for wrongly specified waveform

phi_s = @(t) exp(-abs(t)/s/2) .* sin(3 * 2*pi*t/T).^2;
eval (sprintf ("pdf_s = %s;", func2str (pdf)));
eval (sprintf ("hatphi = %s;", func2str (hatphi)));

hphio = hatphi (to, phi_s,pdf_s);
hphi  = hatphi (t, phi_s,pdf_s);

n_ls = yo.' * hphio ./ sumsq (hphio);
y_ls = n_ls * hphi;
[n_gp y_gp] = gp_estim (phi_s, pdf_s, t, to, io, yo, hphi, hphio);

plot_result (5, t, to, y, yo,y_ls,y_gp,nd,n_ls,n_gp)
** MSE Error **
LS:1.4	53/30
GP:0.48	50/30
s_delayed_decomp-4.png