lans_gammaicdf1.m

模式识别工具包

%	lans_gammaicdf1	- Inverse Gamma cdf (1-D)
%
%	[x]	= lans_gammaicdf1(p,alpha,beta)
%
%	_____OUTPUT_____________________________________________________________
%	x	x value giving cdf p				(vector)
%
%	_____INPUT______________________________________________________________
%	p	percentile	0<p<1				(vector)
%	alpha	value	alpha>-1				(scalar)
%	beta	value	beta>0					(scalar)
%
%	_____EXAMPLE____________________________________________________________
%	lans_gammaicdf1(1,1,1)=
%
%	_____NOTES______________________________________________________________
%	modified from gaminv.m in the MATLAB statistic toolbox
%
%	_____SEE ALSO___________________________________________________________
%	lans_gammacdf1	lans_gammapdf1
%
%	(C) 2000.01.21 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

function [x]	= lans_gammaicdf1(p,alpha,beta);

% Newton's Method

maxiter	= 100;
iter	= 0;

x	= zeros(size(p));
wzero	= find(p==0);
x(wzero)= zeros(size(wzero));
wones	= find(p==1);
x(wones)= ones(size(wones))*inf;
fullx	= x;
therest	= find(p~=0&p~=1);
prest	= p(therest);

if ~any(prest)
	return;
end
	
%	initial guess

mn	= (alpha+1)*beta;
v	= mn*beta;
temp	= log(v+mn*mn);
mu	= 2 * log(mn) - 0.5 * temp;
sigma	= -2 * log(mn) + temp;

x	= exp(lans_gaussicdf1(prest,mu,sigma));

h	= ones(size(prest)); 

% Break out of the iteration loop for three reasons:
%  1) the last update is very small (compared to x)
%  2) the last update is very small (compared to sqrt(eps))
%  3) There are more than 100 iterations. This should NEVER happen. 

while(any(abs(h) > sqrt(eps)*abs(x))  &  max(abs(h)) > sqrt(eps)    ...
                                 & iter < maxiter), 
	iter	= iter + 1;
	h	= (lans_gammacdf1(x,alpha,beta)-prest)./lans_gammapdf1(x,alpha,beta);
	xnew	= x - h;

% Make sure that the current guess stays greater than zero.
% When Newton's Method suggests steps that lead to negative guesses
% take a step 9/10ths of the way to zero:
	ksmall = find(xnew < 0);

	if any(ksmall),
		xnew(ksmall) = x(ksmall) / 10;
		h	= x-xnew;
	end
	x	 = xnew;
end

% Store the converged value in the correct place

if iter == maxiter, 
    fprintf('\nWarning: GAMINV did not converge.\n');
    str = 'The last step was:  ';
    outstr = sprintf([str,'%13.8f'],h);
    fprintf(outstr);
end

xin	= x;
x	= fullx;
x(therest)=xin;