MCMC - simple experiment


In this simple experiment we use a harmonic oscillator ssm model to simulate date from an input and then we use MCMC to recover the parameters of the model.

Contents

Create data


  fs    = 10;
  nsecs = 300;
  fin = 0.2

  % Create a sinusoidal input and arrange it to matrix object
  in = ao(plist('waveform', 'sine wave', 'fs', fs, 'nsecs', nsecs, 'f', fin);  
  min = matrix(in,plist('shape',[1 1]));
 % we do the same for the noise.  
  ns = ao(plist('waveform', 'noise', 'sigma', 0.1, 'nsecs', nsecs, 'fs', fs));  
  mnoise = matrix(ns,plist('shape',[1 1]))

Simulate the data


% declare some variables here
damp = 0.1;
k = 0.1;
values = [k , damp];

% define inNames  and outNames for ssm model. Also declare parameters to be fitted.
inNames = {'COMMAND.force'};
outNames = {'HARMONIC_OSC_1D.position'};
params = {'DAMP','K'}; 

% The model to simulate the inputs specified above.
mod = ssm(plist('built-in','HARMONIC_OSC_1D',...
    'Version','Fitting',...
    'Continuous',1,...
    'SYMBOLIC PARAMS',params));

% Set the parameters k and damp to the ssm model.
mod.setParameters(plist('names',params,'values',values));

% Modify time step and make the model numerical.
mod.modifyTimeStep('newtimestep', 0.1);  

% Generate covariance
mod.generateCovariance;

% Simulate the output after defining the plist.  
pl = plist('Nsamples', fs*nsecs, ...
          'aos variable names',inNames,...
          'aos', min.getObjectAtIndex(1),...
          'return outputs', outNames,...
          'cpsd variable names', cov.find('names'), ...
          'cpsd', cov.find('cov'), ...
          'displayTime', true);
        
[out, plout] = mod.simulate(pl);  
  
  

Declare our model


model = ssm(plist('built-in','HARMONIC_OSC_1D',...
    'Version','Fitting',...
    'Continuous',1,...
    'SYMBOLIC PARAMS',params));

% if we type model in the terminal we will see 
% all the information available for the model we just created

M:   running display
------ ssm/1 -------
      amats: {  [2 x2 ]  }  [1x1]
      bmats: {  [2 x1 ]    []    }  [1x2]
      cmats: {  [1 x2 ]  }  [1x1]
      dmats: { [0] [1] }  [1x2]
   timestep: 0
     inputs:  [1x2 ssmblock]
         1 : COMMAND | force [kg m s^(-2)]
         2 : NOISE | readout [m]
     states:  [1x1 ssmblock]
         1 : HARMONIC_OSC_1D | x [m], xdot [m s^(-1)]
    outputs:  [1x1 ssmblock]
         1 : HARMONIC_OSC_1D | position [m]
  numparams: (M=1)
     params: (K=1, DAMP=1)
    Ninputs: 2
 inputsizes: [1 1]
   Noutputs: 1
outputsizes: 1
    Nstates: 1
 statesizes: 2
 Nnumparams: 1
    Nparams: 2
isnumerical: false
       hist: ssm.hist 
   procinfo: []
   plotinfo: []
       name: HARMONIC_OSC_1D
description: Harmonic oscillator
    mdlfile: 
       UUID: b94172bc-9851-42cd-b56a-ca05503165ad
--------------------

Calculate covariance


ranges = [1e-8 1e-8 ;
          0.5  0.5];
      
pl = plist('FitParams',params,...
  'paramsValues',values,...   
  'inNames',inNames,...
  'ngrid',20,...
  'stepRanges',ranges,...
  'outNames',outNames,...
  'noise',mnoise,...
  'f1',1e-4,...
  'f2',1,...
  'model',mod,...
  'Navs',5,...
  'pinv',false);

mcrb = crb(min,pl);

Fit with MCMC


% First define some variables.
c = 1;
ranges = {[-3 3] [-3 3]};
h = 3;
Tc = [20 80];
numsampl = 200;

% Then we create the plist for the mcmc.
pl = plist('N',numsampl,...
  'J',1,...
  'cov',single(c*mcrb.y),... 
  'range',ranges,...
  'Fitparams',params,...
  'input',min,...
  'noise',mnoise,...
  'model',model,...
  'inNames',inNames,...
  'outNames',outNames,...
  'frequencies',[1e-4 0.5],...
  'fsout',0.2,...
  'Navs',5,...
  'search',true,...
  'Tc',Tc,'heat',h,...
  'jumps',[2e0 1e1 5e2 1e3],...
  'x0',values,...
  'simplex',false,...
  'plot',[1 2],...
  'SNR0',54,...
  'debug',true);

[b smplr]= mcmc(out,pl);



LTP Team