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 parametersm_in = a.e;v_in = a.var;% The order of approximation and related abscissas and weightsn = 3;xi = [0, sqrt(6)/2, -sqrt(6)/2];wi = [2/3, 1/6, 1/6];% Basis points with input var as variancev_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 varianceev_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 pointsf_varvals = feval(nonlin, v_args);% No need to evaluate at relative origin twicef_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 meanm_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;