Home > wafo > wstats > wtcdf.m

wtcdf

PURPOSE ^

Student's T cumulative distribution function

SYNOPSIS ^

F = wtcdf(x,df)

DESCRIPTION ^

 WTCDF  Student's T  cumulative distribution function 
  
  CALL:  F = wtcdf(x,df); 
  
     F = distribution function evaluated at x 
     x = matrix 
    df = degrees of freedom (1,2,....) 
  
  Example: 
    x = linspace(-5,5,200); 
    p1 = wtcdf(x,1); p2 = wtcdf(x,5); 
    plot(x,p1,x,p2)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function F = wtcdf(x,df) 
002 %WTCDF  Student's T  cumulative distribution function 
003 % 
004 % CALL:  F = wtcdf(x,df); 
005 % 
006 %    F = distribution function evaluated at x 
007 %    x = matrix 
008 %   df = degrees of freedom (1,2,....) 
009 % 
010 % Example: 
011 %   x = linspace(-5,5,200); 
012 %   p1 = wtcdf(x,1); p2 = wtcdf(x,5); 
013 %   plot(x,p1,x,p2) 
014  
015 % tested on matlab 5.3 
016 %History: 
017 %revised pab 22.05.2003 
018 % -added new methods for df==1 or df==2 and for region1= x^2<df^2 
019 %revised pab 29.10.2000 
020 % adapted from stixbox 
021 % -added nargchk, comnsize, mxdf +  check on floor(df)==df 
022 %by      Anders Holtsberg, 18-11-93 
023 %       Copyright (c) Anders Holtsberg 
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 %df = min(df,1000000); % make it converge and also accept Inf. 
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

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