update_params.m

Nonlinear dynamical factor analysis Matlab package

function params = update_params(fs, tfs, sources, data, net, tnet, params),
% UPDATE_PARAMS  Update estimates for different parameters

% Copyright (C) 2002-2005 Harri Valpola, Antti Honkela and Matti Tornio.
% 
% This package comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package.  This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.

params.noise = estimatevars(fs - data, params.hyper.noise, params.noise);
tmpsrc = probdist(sources.e - tfs.e, sources.nvar + tfs.var);
params.src = estimatevars(tmpsrc, params.hyper.src, params.src);
params.net.w2var = estimatevars(net.w2, params.hyper.net.w2var, ...
				params.net.w2var, 1);
params.tnet.w2var = estimatevars(tnet.w2, params.hyper.tnet.w2var, ...
				    params.tnet.w2var, 1);
[params.hyper.net.w2var.mean, params.hyper.net.w2var.var] = ...
    estimatemeanvars(params.net.w2var, params.prior.net.w2var.mean, ...
		     params.prior.net.w2var.var, ...
		     params.hyper.net.w2var.var);
[params.hyper.tnet.w2var.mean, params.hyper.tnet.w2var.var] = ...
    estimatemeanvars(params.tnet.w2var, params.prior.tnet.w2var.mean, ...
		     params.prior.tnet.w2var.var, ...
		     params.hyper.tnet.w2var.var);
[params.hyper.noise.mean, params.hyper.noise.var] = ...
	estimatemeanvars(params.noise, params.prior.noise.mean, ...
			 params.prior.noise.var, params.hyper.noise.var, 1);
[params.hyper.net.b1.mean, params.hyper.net.b1.var] = ...
    estimatemeanvars(net.b1, params.prior.net.b1.mean, ...
		     params.prior.net.b1.var, params.hyper.net.b1.var, 1);
[params.hyper.net.b2.mean, params.hyper.net.b2.var] = ...
    estimatemeanvars(net.b2, params.prior.net.b2.mean, ...
		       params.prior.net.b2.var, params.hyper.net.b2.var, 1);
[params.hyper.tnet.b1.mean, params.hyper.tnet.b1.var] = ...
    estimatemeanvars(net.b1, params.prior.tnet.b1.mean, ...
		     params.prior.tnet.b1.var, params.hyper.tnet.b1.var, 1);
[params.hyper.tnet.b2.mean, params.hyper.tnet.b2.var] = ...
    estimatemeanvars(net.b2, params.prior.tnet.b2.mean, ...
		     params.prior.tnet.b2.var, params.hyper.tnet.b2.var, 1);
[params.hyper.src.mean, params.hyper.src.var] = ...
    estimatemeanvars(params.src, params.prior.src.mean, ...
		     params.prior.src.var, params.hyper.src.var, 1);



function v_x = estimatevars(x, vprior, v_x0, dim, missing)
% ESTIMATEVARS  Estimate parameters for variables with zero mean
%
%    ESTIMATEVARS can be used to estimate variance parameter
%    values for variables that are assumed to have zero mean.
%    V_X = ESTIMATEVARS(X, VPRIOR, V_X0, DIM) finds the estimate V_X for
%    variance parameters of X with previous value V_X0 and
%    hyperparameter value VPRIOR.  Different samples of data
%    are assumed to be along dimension DIM (default: 2).

% Copyright (C) 2002-2004 Harri Valpola, Antti Honkela and Matti Tornio.
%
% This package comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package.  This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.

if nargin < 4
  dim = 2;
end

epsilon = 1e-5;
minstep = -0.5;
basex = v_x0.e;

N = size(x, dim);

% Discard unknown data values
if nargin >= 5,
  x(find(missing)) = probdist(0, 0);
  N = N - sum(missing);
end

sueff = exp(2*vprior.var.e - 2*vprior.var.var);
xval = sum(x.e .^ 2 + x.var, dim);

beta = sueff .* xval .* exp(-2 * basex + 2 * v_x0.var);
gamma = vprior.mean.e - N * sueff - basex;

t = zeros(size(v_x0));
% solve t - beta * exp(-2 t) - gamma = 0
% using Newton's iteration

step = ones(size(v_x0)) + epsilon;

while max(abs(step)) > epsilon
  step = (t - beta .* exp(-2 * t) - gamma) ./ (-1 - 2*beta.*exp(-2 * t));
  step = step .* (step >= minstep) + minstep * (step < minstep);
  t = t + step;
end

new_mean = basex + t;
new_var = sueff ./ (1 + 2*(t - gamma));

v_x = probdist(new_mean, new_var);



function [m_x, v_x] = estimatemeanvars(x, mprior, vprior, v_x0, dim)
% ESTIMATEMEANVARS  Estimate parameters for variables with nonzero
%                   mean and variance
%
%    ESTIMATEMEANVARS can be used to estimate parameter values
%    for variables that are assumed to have nonzero mean and variance.

%    [M_X, V_X] = ESTIMATEMEANVARS(X, MPRIOR, VPRIOR, V_X0, DIM) finds
%    the estimates M_X and V_X for parameters of X with previous value
%    of V_X V_X0 and priors for the parameters MPRIOR and VPRIOR.
%    Different samples of data are assumed to be along dimension DIM
%    (default: 2).

% Copyright (C) 2002 Harri Valpola and Antti Honkela.
%
% This package comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package.  This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.

if nargin < 5
  dim = 2;
end

N = size(x, dim);
veff = exp(2*v_x0.e - 2*v_x0.var);
vmeff = exp(2*mprior.var.e - 2*mprior.var.var);
new_var = 1./(N./veff + 1/vmeff);
new_mean = (sum(x.e, dim)./veff + mprior.mean.e/vmeff).*new_var;
m_x = probdist(new_mean, new_var);

epsilon = 1e-5;
minstep = -0.5;
basex = v_x0.e;

vveff = exp(2*vprior.var.e - 2*vprior.var.var);
xval = sum(x.e .^ 2 + x.var, dim) + N*(new_var + new_mean.^2) - ...
    2*sum(x.e, dim).*new_mean;

beta = vveff .* xval .* exp(-2 * basex + 2 * v_x0.var);
gamma = vprior.mean.e - N * vveff - basex;

t = zeros(size(v_x0));
% solve t - beta * exp(-2 t) - gamma = 0
% using Newton's iteration

step = ones(size(v_x0)) + epsilon;

while max(abs(step)) > epsilon
  step = (t - beta .* exp(-2 * t) - gamma) ./ (-1 - 2*beta.*exp(-2 * t));
  step = step .* (step >= minstep) + minstep * (step < minstep);
  t = t + step;
end

new_mean = basex + t;
new_var = vveff ./ (1 + 2*(t - gamma));

v_x = probdist(new_mean, new_var);