Home > wafo > wstats > wbetainv.m

wbetainv

PURPOSE ^

Inverse of the Beta distribution function

SYNOPSIS ^

x = wbetainv(F,a,b)

DESCRIPTION ^

 WBETAINV  Inverse of the Beta distribution function
 
  CALL:  x = wbetainv(F,a,b)
 
   x   = inverse cdf for the Beta distribution evaluated at F.
  a, b = distribution parameters
 
  Example:
    F = linspace(0,1,100);
    x = wbetainv(F,1,2);
    plot(F,x)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function x = wbetainv(F,a,b)
002 %WBETAINV  Inverse of the Beta distribution function
003 %
004 % CALL:  x = wbetainv(F,a,b)
005 %
006 %  x   = inverse cdf for the Beta distribution evaluated at F.
007 % a, b = distribution parameters
008 %
009 % Example:
010 %   F = linspace(0,1,100);
011 %   x = wbetainv(F,1,2);
012 %   plot(F,x)
013 
014 % tested on matlab 5.3
015 %History:
016 %revised pab 29.10.2000
017 % adapted from stixbox
018 % -added nargchk, comnsize
019 % -refined the Newton-Raphson method
020 %       Anders Holtsberg, 27-07-95
021 %       Copyright (c) Anders Holtsberg
022 
023 error(nargchk(3,3,nargin))
024 [errorcode F,a,b] = comnsize(F,a,b);
025 if errorcode>0,
026   error('x, a and b must be of common size or scalar');
027 end
028 
029 x = zeros(size(F));
030 
031 ok = (a>0 & b>0);
032 
033 k = find(F>0 & F<1 & ok);
034 if any(k)
035  
036   bk = b(k); %min(b(k),100000);
037   if any(bk>10^5), warning('b is too large'),end
038   ak = a(k);
039   Fk = F(k);
040   xk = max(min(ak ./ (ak+bk),1-sqrt(eps)),sqrt(eps));
041  
042   dx = ones(size(xk));
043   max_count=100;
044   ix = find(xk); iy=0;
045   while (any(ix) & iy<max_count),
046     iy=iy+1;
047     xi=xk(ix);ai=ak(ix);bi=bk(ix);
048     dx(ix) = (wbetacdf(xi,ai,bi) - Fk(ix)) ./wbetapdf(xi,ai,bi);
049     dx(ix)= min(abs(dx(ix)),.5/iy).*sign(dx(ix));
050     xi = xi - dx(ix);
051     
052     % Make sure that the current guess is larger than zero and less than 1
053     xk(ix) = xi + 0.1*(dx(ix) - 9*xi).*(xi<=0) +...
054     0.38*(dx(ix)-6.2*xi +6.2).*(xi>=1) ;
055     
056     ix=find((abs(dx) > sqrt(eps)*abs(xk))  &  abs(dx) > sqrt(eps));
057         disp(['Iteration ',num2str(iy),...
058         '  Number of points left:  ' num2str(length(ix)) ]),
059   end
060  
061   x(k)=xk;
062   if iy == max_count, 
063     warning('wbetainv did not converge.');
064     str = ['The last step was less than:  ' num2str( max(abs(dx(ix))) )];
065     disp(str);
066     k(ix)
067   end
068 end
069 
070 
071 k2=find(F==1&ok);
072 if any(k2)
073   x(k2)=ones(size(k2));
074 end
075 
076 k3=find(~ok);
077 if any(k3)
078   tmp=NaN;
079   x(k3)=tmp(ones(size(k3)));
080 end
081 
082 
083 
084 
085 
086 
087

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003