qinv.m

来自「这是在网上下的一个东东」· M 代码 · 共 30 行

%
%  qinv=qinv(y,a,u,beta)
%
%  This function computes the inverse of the posterior CDF.  It us used
%  by gensol.
%
function qinv=qinv(y,a,u,beta)

% if a = 0, use the limit

if (abs(a) < eps*10^8)
	qinv=u*y;

% if beta << -1, use the asymptotic approximation
elseif (beta < -300)

	term1=log(y*a^2*(3*beta-a)/(-6*beta^3)+exp(-a/beta));
	mnew=u;
	stopcrit=1;
	while stopcrit > eps*10^8
		mold=mnew;
		fofx=log(a*(mold+1/a-u-1/beta))+a*(u-mold)-term1;
		fprime=-a+1/(mold+1/a-u-1/beta);
		mnew=mold-fofx/fprime;
		stopcrit=abs((mold-mnew)/mold);
	end
	qinv=mnew;
else
	qinv=-log(1-y+y*exp(-a*u))/a;
end