001 function x = wbetainv(F,a,b)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
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);
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
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