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:
 comnsize Check if all input arguments are either scalar or of common size. wbetacdf Beta cumulative distribution function wbetapdf Beta probability density function error Display message and abort function. nan Not-a-Number. num2str Convert number to string. (Fast version) warning Display warning message; disable or enable warning messages.
This function is called by:
 wbetarnd Random matrices from a Beta distribution wfinv Inverse of the Snedecor's F distribution function wtinv Inverse of the Student's T distribution function

## 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
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