nonlin_hermite.m

非线型因素分析matlab仿真程序包

function [b, aux] = nonlin_hermite(a, nonlin)
% TANH Hyperbolic tangent of a propability distribution
%
% The result is calculated using Gauss-Hermite quadrature formula of
% order 3 for other parts, except for multivars which are handled by
% linearising the nonlinearity at mean
%

% Copyright (C) 1999-2004 Antti Honkela, Harri Valpola,
% and Xavier Giannakopoulos.
%
% 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.

% The input parameters
m_in = a.e;
v_in = a.var;

% The order of approximation and related abscissas and weights
n = 3;
xi = [0, sqrt(6)/2, -sqrt(6)/2];
wi = [2/3, 1/6, 1/6];

% Basis points with input var as variance
v_args = zeros([size(m_in), n]);
argtemp = sqrt(2 * v_in);
for k=1:length(xi),
  v_args(:, :, k) = m_in + xi(k) * argtemp;
end
%v_args = gplus(m_in, gtimes(reshape(xi, [1, 1, n]), argtemp));

% Basis points with extravar as variance
ev_args = zeros([size(m_in), n]);
argtemp = sqrt(2 * a.extra);
for k=1:length(xi),
  ev_args(:, :, k) = m_in + xi(k) * argtemp;
end
%ev_args = gplus(m_in, gtimes(reshape(xi, [1, 1, n]), argtemp));
clear argtemp

% Nonlinearity evaluated at required points
f_varvals = feval(nonlin, v_args);

% No need to evaluate at relative origin twice
f_evvals = f_varvals;
f_evvals(:, :, 2:end) = feval(nonlin, ev_args(:, :, 2:end));

% Components of the sum to evaluate output mean with var as
% input variance
%v_sumparts = repmat(reshape(wi, [1, 1, n]), [size(m_in), 1]) .* f_varvals;
v_sumparts = gtimes(reshape(wi, [1, 1, n]), f_varvals);

% Output mean
m_out = sum(v_sumparts, 3);

% Function deviations
%f_vardevs = f_varvals - repmat(m_out, [1, 1, n]);
%f_evdevs = f_evvals - repmat(m_out, [1, 1, n]);
f_vardevs = gminus(f_varvals, m_out);
f_evdevs = gminus(f_evvals, m_out);

% Output var and extravar
%v_out = sum(repmat(reshape(wi, [1, 1, n]), [size(m_in), 1]) .* ...
%	    (f_vardevs .^ 2), 3);
%ev_out = sum(repmat(reshape(wi, [1, 1, n]), [size(m_in), 1]) .* ...
%	     (f_evdevs .^ 2), 3);
v_out = sum(gtimes(reshape(wi, [1, 1, n]), f_vardevs .^ 2), 3);
ev_out = sum(gtimes(reshape(wi, [1, 1, n]), f_evdevs .^ 2), 3);

% Components of the sum to evaluate output mean with extravar as
% input variance
%ev_sumparts = repmat(reshape(wi, [1, 1, n]), [size(m_in), 1]) .* ...
%	 feval(nonlin, ev_args);

% Normalised sum components
%v_sumparts0 = ...
%    v_sumparts - repmat(reshape(wi, [1, 1, n]), [size(m_out), 1]) .* ...
%    repmat(m_out, [1, 1, n]);
%ev_sumparts0 = ...
%    ev_sumparts - repmat(reshape(wi, [1, 1, n]), [size(m_out), 1]) .* ...
%    repmat(m_out, [1, 1, n]);

d = sqrt(v_out ./ v_in);

% Multiply the multivar vectors, i.e. calculate
%  m(:,i,:) = d .* a.multivar(:,i,:);
[dim1 dim2 dim3] = size(a.multi);
%new_mv = repmat(reshape(d, [dim1 1 dim3]), [1 dim2 1]) .* a.multivar;
new_mv = gtimes(reshape(d, [dim1 1 dim3]), a.multi);

b.e = m_out;
b.var = v_out;
b.multi = new_mv;
b.extra = ev_out;

aux.f_varvals = f_varvals;
aux.f_vardevs = f_vardevs;
aux.f_evvals = f_evvals;
aux.f_evdevs = f_evdevs;
aux.v_sumparts = v_sumparts;
aux.v_args = v_args;
aux.ev_args = ev_args;