getfk.m

这是在网上下的一个东东

%%%%%%%%%% start of function getfk.m %%%%%%%%%%
%
% Function to calculate the F vector (Equation 2.28)
%
% Written by: Roseanna M. Neupauer
% Modification Date: April 24, 1999
% Modification Date: June 25, 1999 (added option to use additive noise,
%		multiplicative noise, or a combination of the two)
% Modification Date: March 31, 2000 (corrected the asymptotic appoximation
%		for beta < -blarge)
%
% [fe,mhat,a]=getfk(lambda,g,beta,upper,d,noise,...
%          nearzero,large)
%
% Inputs
%     lambda       array containing the Lagrange multipliers
%     g            matrix of scaled kernel functions
%     beta         array of Lagrange multipliers
%     upper        array containing upper limit of prior 
%                    distributions
%     d            array containing observed data.
%     noise        noise level
%     nearzero     value below which the asymptotic 
%                    approximation to zero is used
%     large        value above which the asymptotic 
%                    approximation to infinity is used
%
% Outputs
%     fe           array containing the F vector
%     mhat         array containing the expected value of the
%                    model, based on the posterior distribution
%     a            array containing the vector, a

function [fe,mhat,a]=getfk(lambda,g,beta,upper,d,...
     noise,nearzero,large)

% define array sizes and create new arrays
[ndata,nt]=size(g);
fe=zeros(ndata,1);
b=zeros(nt,1);
a=zeros(nt,1);
mhat=zeros(nt,1);
blarge=100.;

% calculate the a vector (page 35)
for k=1:nt
     aa=0.;
     for l=1:ndata
          aa=aa+lambda(l)*g(l,k);
     end
     a(k)=beta(k)+aa;
     b(k)=aa;
end

% calculate mhat 
mau=a.*upper;
for j=1:nt 
     if mau(j) == 0.

          % Equation 2.28          
          mhat(j)=upper(j)/2;

     elseif (beta(j) < -blarge)

          % Equation 3.11
%          mhat(j)=2*upper(j)/3-1/(3*beta(j))*(1-upper(j)*...
%               (beta(j)+b(j)))+...
%               2/(9*beta(j)^2)*(-beta(j)-b(j)+upper(j)*b(j)*...
%               (beta(j)+b(j)/2));
	mhat(j)=upper(j)+(b(j)-beta(j))/beta(j)/(b(j)-3*beta(j));
	  a(j)=b(j);

     elseif (abs(mau(j)) < nearzero) 
          au=a(j)*upper(j);

          % Equation 3.12
          mhat(j)=upper(j)*(12.-8.*au+3.*au^2)/...
               (24.-12.*au+4.*au^2-au^3);

     elseif mau(j) < -1*10^2

          % page 60
          mhat(j)=upper(j)+1/a(j);

     elseif (mau(j) > 10^4)

          % page 60
          mhat(j)=1/a(j);

     else

          % Equation 2.28
          mhat(j)=(exp(-mau(j)).*(-mau(j)-1)+1)./...
               (a(j).*(1-exp(-mau(j))));

     end     
end

% calculate the F matrix 
errfactor=(sqrt(ndata)*noise(1)+norm(d)*noise(2))/norm(lambda);

% Equations 3.10 and B.11
for k=1:ndata
     ff=d(k);
     for l=1:nt
          ff=ff-g(k,l)*mhat(l);
     end
     fe(k)=ff+errfactor*lambda(k);
end

return
%%%%%%%%%%  end of function getfk.m  %%%%%%%%%%