001 function F = wtcdf(x,df)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026 error(nargchk(2,2,nargin))
027 [errorcode x,df] = comnsize(x,df);
028 if errorcode>0,
029 error('x and df must be of common size or scalar');
030 end
031
032 F = zeros(size(x));
033
034
035 mxdf = 10^6;
036
037 k1 = find(df==1 );
038 if any(k1)
039 F(k1) = ( 1 + 2*atan(x(k1))/pi )/2;
040 end
041 k2 = find(df==2);
042 if any(k2)
043 x2= x(k2);
044 F(k2) = ( 1 + x2./sqrt( 2 + x2.*x2 ))/2;
045 end
046
047 ok = (0<df & df==floor(df));
048 region1 = (abs(x)<sqrt(abs(df)));
049 k3 = find(ok & region1 & (2 < df) & (df<mxdf) );
050 if (any(k3)),
051 dfk = df(k3);
052 xk = x(k3);
053 xk2 = xk.*xk;
054 cssthe = 1./( 1 + xk2./dfk );
055 nuVec = unique(dfk(:));
056 Fk = zeros(size(dfk));
057 for nu = nuVec(:).'
058 knu = find(dfk==nu);
059 Fk(knu) = evalPoly(nu,xk(knu),xk2(knu),cssthe(knu));
060 end
061 F(k3) = Fk;
062 end
063
064 k = find(ok & ~region1 & (2<df) & df<mxdf );
065 if any(k),
066 neg = x(k)<0;
067 tmp = 1-(1-wfcdf(x(k).^2,1,df(k)))/2;
068 F(k) = tmp + (1-2*tmp).*neg;
069 end
070
071 k1=find(ok & df>=mxdf);
072 if any(k1)
073 F(k1) = wnormcdf(x(k1),0,1);
074 end
075
076
077 k2 = find(~ok);
078 if any(k2)
079 F(k2)=NaN;
080 end
081
082 return
083 function F = evalPoly(nu,t,tt,cssthe)
084
085 polyn = 1;
086 for j = nu-2 : -2 : 2
087 polyn = 1 + ( j - 1 )*cssthe.*polyn/j;
088 end
089 if ( mod( nu, 2 ) == 1 )
090 ts = t/sqrt(nu);
091 F = ( 1 + 2*( atan(ts) + ts.*cssthe.*polyn )/pi )/2;
092 else
093 snthe = t./sqrt( nu + tt );
094 F = ( 1 + snthe.*polyn )/2;
095 end
096 F = max(0, min(F, 1) );
097 return