vssminpn.m

基于变分贝叶斯方法的状态空间模型学习程序

%|              Variational Bayesian State-Space models                 |
%|                                                                      |
%|      Copyright Matthew J. Beal 22/05/01 version 3.0 (11/08/03)       |
%|              http://www.cs.toronto.edu/~beal/software                |
%
%USAGE:  net = vssminp(Yn,inpn,k,its,dispopt,cFbool,hypit,net);
%
%-Reqrd-
% Yn      N obs sequences    cell (1 by N), each cell is matrix (p by Tn)
% inpn    input sequence(s)  as for Yn, but can also be just one sequence
%
%-Optnl-  description (df default setting)              (possible values)
% k       maximum hidden state dimension (df p)              (at least 1)
% its     maximum it#,iterations (df 100)                    (at least 1)
% dispopt verbosity of output (0=nothing, df 1)                 (0,1,2,3)
% cFbool  calc F @each it (df 0, but F always calc'd at last it)    (0,1)
% hypit   it# to optimise dyn,oput,&other hyps (df [10 20 30])   (1 by 3)
% net     a previous output (overwrites k setting)     (matlab structure)
%
%Returns 'net', a structure consisting of:
%
% Optim. hyperparams: net.hyp.{pa,pb,alpha,beta,gamma,delta,X0m_p,X0ci_p}
% Sufficient statistics of the parameter posterior:  net.param.{SigA,...}
% Expected natural parameters:                    net.exp.{A,B,AA,AB,...}
% Sufficient stastistics of the hidden state:
%                             net.hidden.{X0mn,Xmn,X0cin,Xcin,Ups0n,Upsn}
% Histories/trajectories of the lower bound and hyperparameters:
%                                            net.hist.{F,pa,pb,alpha,...}
% E.g. net.exp.B is the mean dynamics matrix, and
%      net.param.SigB is the covariance of each of its rows.
% Note: not the same for matrix A! (see ch.5 of my thesis)
%
%v3.0 Matthew J. Beal GCNU 01/03/03

function net = vssminpn(Yn,inpn,k,its,dispopt,cFbool,hypit,net);

% Initialise dimensions of inputs/outputs
n = size(Yn,2);			% Yn is a cell array of size ( 1 x #sequences )
Tn = zeros(1,n);		% Tn is a row vector containing length of each seq
for i=1:n
  Tn(i) = size(Yn{i},2);	% sequence length may vary
end
p = size(Yn{1},1);		% take first obs seq to set the observ. dim
pinp = size(inpn{1},1);		% take first input sequence to set input dim

% If just one input sequence is specified, replicate it to be used as the
% input sequence for every sequence.
if size(inpn,2)~=n		% then use the first input sequence as master
  mstseq=inpn{1};
  if size(mstseq,2)<max(Tn)
    fprintf('Error: I will only replicate the first inp seq, and it is not long enough'); return;
  end
  fprintf('[ replicating first inp seq ]\n');
  for i=1:n
    inpn{i}=mstseq(:,1:Tn(i));
  end
end

% Compute some statistics of the data & inputs (no need to put in VBEM loop)
Ydd = zeros(p,p); Udd=zeros(pinp,pinp); UY=zeros(pinp,p);
for i=1:n
  Ydd = Ydd+Yn{i}*Yn{i}';
  Udd = Udd+inpn{i}*inpn{i}';
  UY = UY+inpn{i}*Yn{i}';
end

if nargin<3, k = p; end		% if k not specified, set equal to dimensions of output
				%   this irrelevant if 'net' is provided in function call
if nargin<4, its = 100; end	% if no limit specified then set to 100 iterations
if nargin<5, dispopt=2; end	% degree of verbosity set to just print lower bound progress (no plots)
if nargin<6, cFbool=0; end      % default is to leave F uncalculated (faster)
if nargin<7,
  dynhypit=10;			% it# after which dynamics hyps optimised
  outhypit=20;			% it# after which output hyps optimised
  etchypit=30;			% it# after which other (input & noise hyps) are optimised
 else
  dynhypit=hypit(1); outhypit=hypit(2); etchypit=hypit(3); % user-specified schedule
end

if nargin<8,			% Set up all hyperparameters and initialisation right now
 
  if dispopt>=1; fprintf('Initialising hyps,'); end
  X0m_p = zeros(k,1);		% time zero mean of X set to zero
  X0ci_p = eye(k,k);		% time zero precision X set to identity
  pa = 1; pb = 1;		% shape and precision of noise precision prior
  alpha = 1*ones(k,1);		% precision vector, each element for a COL of A
  beta = 1*ones(pinp,1);	% precision vector, each element for a COL of B
  gamma = 1*ones(k,1);		% pseudo-precision vector, each element for a COL of C
  delta = 1*ones(pinp,1);	% pseudo-precision vector, each element for a COL of D

  % Initialise a randomised hidden state
  %   This is somehow preferable to initialising a set of expected natural parameters
  %   from their prior, which technically speaking is the more correct thing to do --- Sorry!
  if dispopt>=1; fprintf(' and randomly instantiating hidden state sequences'); end
  for i=1:n
    X0mn{i} = X0m_p;				% hidden state means
    Xmn{i} = randn(k,Tn(i));
    X0cin{i} = X0ci_p;				% hidden state inverse covariances
    Xcin{i} = repmat(eye(k,k),[1 1 Tn(i)]);
    Ups0n{i} = eye(k,k);			% hidden state cross-correlations
    Upsn{i} = repmat(eye(k,k),[1 1 Tn(i)]);
  end
  % Set up hyperparameter histories
  pahist = pa; pbhist = pb; alphahist = alpha; betahist = beta; gammahist = gamma; deltahist = delta; X0m_phist = X0m_p; X0ci_phist = X0ci_p;
  hist = -Inf*ones(1,7);

else				% Use the net provided in the function call

  if dispopt>=1; fprintf('Extracting hyperparameters, their histories, and parameter posterior  from previously learnt network'); end
  % Extract the hyperparameters...
  pa	= net.hyp.pa;
  pb	= net.hyp.pb;
  alpha	= net.hyp.alpha; k=size(net.hyp.alpha,1);
  beta	= net.hyp.beta;
  gamma	= net.hyp.gamma;
  delta	= net.hyp.delta;
  X0m_p	= net.hyp.X0m_p;
  X0ci_p= net.hyp.X0ci_p;
  % Extract also the expected natural parameters...
  A	= net.exp.A;
  B	= net.exp.B;
  AA	= net.exp.AA;
  AB	= net.exp.AB;
  BB	= net.exp.BB;
  rho	= net.exp.rho;
  lnrho	= net.exp.lnrho;
  CrhoC	= net.exp.CrhoC;
  rhoC	= net.exp.rhoC;
  CrhoD	= net.exp.CrhoD;
  rhoD	= net.exp.rhoD;
  DrhoD	= net.exp.DrhoD;
  % Finally now infer the hidden state sequence in the 0th VBEM step.
  %   NOTA BENE: this way it is possible to use the posterior over parameters 
  %   for training set A, and then continue training FROM THAT POINT on a second 
  %   training set B (which may have sequences of different length).
  % VBE step : q(X) is Gaussian
  for i=1:n
    [lnZpTn{i},Xm_a,Xci_a] = forwardpass(Yn{i},inpn{i},X0ci_p,X0m_p,A,AA,AB,B,BB,rho,lnrho,CrhoC,rhoC,CrhoD,rhoD,DrhoD,cFbool);
    [Xm_b,Xci_b,X0m_b,X0ci_b] = backwardpass(Yn{i},inpn{i},A,AA,B,AB,CrhoC,rhoC,CrhoD);
    [Xmn{i},Xcin{i},Upsn{i},X0mn{i},X0cin{i},Ups0n{i}] = getmarginals(Xm_a,Xci_a,Xm_b,Xci_b,X0m_p,X0ci_p,X0m_b,X0ci_b,A,AA,CrhoC);  
  end  
  % Gather up previous hyperparameter histories
  pahist = net.hist.pa;
  pbhist = net.hist.pb;
  alphahist = net.hist.alpha;
  betahist = net.hist.beta;
  gammahist = net.hist.gamma;
  deltahist = net.hist.delta;
  X0m_phist = net.hist.X0m_p;
  X0ci_phist = net.hist.X0ci_p;
  hist = net.hist.F;

end

% Display some summary information on the data sequences and proposed model size
if dispopt>=1;
  fprintf('.\nSummary:');
  fprintf(' # of sequences           : %i\n',n);
  fprintf('         Min,Max sequence lengths : %i,%i\n',min(Tn),max(Tn));
  fprintf('         Total # of observations  : %i\n',sum(Tn,2));
  fprintf('         # hidden states in model : %i',k);
  if ~cFbool, fprintf('\nSkipping F calculations (fast, but no history)'); end
end

% Set up figure plots:
if dispopt>=3,
  hf = figure; clf
  subplot(211);
    hFpart = plot(zeros(2,6)); title('F constituents');
    legend('B','A|B','rho','D|rho','C|rho,D','F(Y)');
  subplot(212);
    hFdiff = plot(zeros(1,1)); title('log rate of change of F');
    legend('log \DeltaF(Y)');
  hf2 = figure; clf; 
  subplot(231); halpha = plot(zeros(2,k)); title('prior log var A_{|.}');
  subplot(232); hbeta = plot(zeros(2,pinp)); title('prior log var B_{|.}');
  subplot(234); hgamma = plot(zeros(2,k)); title('prior log var C_{|.}');
  subplot(235); hdelta = plot(zeros(2,pinp)); title('prior log var D_{|.}');
  subplot(233); hrhomean = plot(zeros(2,p)); title('prior \rho mean');
  subplot(236); hrhovar = plot(zeros(2,p)); title('prior \rho log var');
  set(gcf,'doublebuffer','on');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN VBEM iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

it = 0; dF=Inf; tic;
while 1 % note there is a loop break after the F-calculation, of the sort: ~( (it<its) & (dF>1e-14) )

  it = it+1; if (it==4) & (dispopt>=1), fprintf('\nApproxTime for %i iterations: %2.1f mins (''.''=50its)',its,toc/3*(its-3)/60); end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % UPDATE A,B
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % VBM step (dynamics model)
  % Compute parameters to represent Q(A,B)
  WA = zeros(k,k);
  SA = zeros(k,k);
  GA = zeros(k,pinp);
  MB = zeros(k,pinp);
  for i=1:n
    WA = WA+inv(X0cin{i})+ X0mn{i}*X0mn{i}'+ Xmn{i}(:,1:(Tn(i)-1))*Xmn{i}(:,1:(Tn(i)-1))';
    for t = 1:(Tn(i)-1);
      WA = WA+inv(Xcin{i}(:,:,t));
    end
    SA = SA+ Ups0n{i}+X0mn{i}*Xmn{i}(:,1)' + sum(Upsn{i}(:,:,1:(Tn(i)-1)),3)+Xmn{i}(:,1:Tn(i)-1)*Xmn{i}(:,2:Tn(i))';
    GA = GA+ X0mn{i}*inpn{i}(:,1)' + Xmn{i}(:,1:Tn(i)-1)*inpn{i}(:,2:Tn(i))';
    MB = MB+ Xmn{i}*inpn{i}';
  end
  
  SigA = inv( WA + diag(alpha) );