mre.m

这是在网上下的一个东东

%
%  [x,lambda,beta,a,p5,p95]=mre(G,d,upper,lower,expvalue,noise)
%
%  This routine performs minimum relative entropy (MRE) inversion for 
%  a discrete linear inverse problem of the form 
%
%                    Gx=d
%
%  with upper and lower bounds on x,
%
%                    l <= x <= u
%
%  and specified prior expected value expvalue.  
%
%  The user also specifies a noise level.  The MRE procedure will
%  ensure that the posterior mean solution satisifies 
%
%              || Gx-d || <= sqrt(n)*noise(1)+norm(d)*noise(2)
%
%  Note that the routine may fail to converge.  This usually means that 
%  the noise tolerance was too small.  
%
%  Inputs:
%         G                the matrix G in Gx=d
%         d                the data vector d in Gx=d
%         upper            upper bounds on x
%         lower            lower bounds on x
%         expvalue         Expected value of x
%         noise            noise level
%                          noise(1)   absolute noise
%                          noise(2)   relative noise
%
%  Outputs:
%         x                Mean value of the posterior pdf
%         lambda           Lagrange multipliers
%         beta             parameters of the posterior pdf 
%         a                parameters of the posterior pdf          
%         p5               5% probability level of the posterior pdf
%         p95              95% probability level of the posterior pdf
%
function [x,lambda,beta,a,p5,p95]=mre(G,d,upper,lower,expvalue,noise)
%
%  First, subtract off the effect of the lower bound.
%
d2=d-G*lower;
expvalue2=expvalue-lower;
upper=upper-lower;
%
%  Setup default values for the tolerances used by domre.
%
leftbegin=-20;
rightbegin=20;
tolbeta=0.0001;
maxiter=30;
tollam=0.01;
lamiter=30;
tolls=0.00001;
nearzero=0.0001;
large=300;
%
% Call domre
%
[x,lambda,beta,a,p5,p95]=domre(G,d2,upper,expvalue2,noise,leftbegin,...
                                      rightbegin,tolbeta,maxiter,tollam,...
                                      lamiter,tolls,nearzero,large);
%
% Adjust the solution for the lower bounds.
%
x=x+lower;
p5=p5+lower;
p95=p95+lower;
%
% And we're done!
%