statistic.cpp Source File

00001 
00002 #ifndef Statistic_SOURCE
00003 #define Statistic_SOURCE
00004 
00005 #include "statistic.h"
00006 #include <cmath>
00007 
00010 namespace PLib {
00011 
00012   void initStatistic(){
00013     MaximumIterations = 100 ; 
00014   }
00015 
00021   template <class T>
00022     T lnOfGamma(T xx)
00023     {
00024       double x,y,tmp,ser;
00025       static double cof[6] = { 76.18009172947146, -86.50532032941677,
00026                                24.01409824083091, -1.231739572450155,
00027                                0.1208650973866179e-2,-0.5395239384953e-5};
00028       int j;
00029       y=x=xx;
00030       tmp=x+5.5 ;
00031       tmp -= (x+0.5)*log(tmp);
00032       ser=1.000000000190015;
00033       for(j=0;j<=5;++j) ser += cof[j]/++y ;
00034       return T(-tmp+log(2.5066282746310005*ser/x));
00035     }
00036  
00037   template <class T>
00038     T factorial(int n)
00039     {
00040       static int ntop=4;
00041       static T a[33] = {1.0,1.0,2.0,6.0,24.0} ;
00042       if(n<0){
00043 #ifdef USE_EXCEPTION       
00044         throw MatrixInputError();
00045 #else
00046         Error error("factorial");
00047         error << "The input value was lower than 0.\n" ;
00048         error.fatal();
00049 #endif
00050       }
00051       if(n>32)
00052         return exp(lnOfGamma(n+1.0));
00053       while(ntop<n){
00054         int j=ntop++;
00055         a[ntop]=a[j]*ntop ;
00056       }
00057       return a[n] ;
00058     }
00059 
00060   template <class T>
00061     T binomialCoefficient(int n, int k)
00062     {
00063       return (T)floor(0.5+exp(lnOfFactorial<T>(k)-lnOfFactorial<T>(n-k)));
00064     }
00065 
00066   template <class T>
00067     T lnOfFactorial(int n)
00068     {
00069       static T a[101] ;
00070       if(n<0){
00071 #ifdef USE_EXCEPTION       
00072         throw MatrixInputError();
00073 #else
00074         Error error("lnOfFactorial");
00075         error << "The input value was lower than 0.\n" ;
00076         error.fatal();
00077 #endif
00078       }
00079       if(n<=1) return 0;
00080       if(n <= 100) 
00081         return a[n] ? a[n] : (a[n] = lnOfGamma(n+1.0)) ;
00082       return lnOfGamma(n+1.0) ; // out of range of table
00083     }
00084   
00085   template <class T>
00086     T beta(T z, T w)
00087     {
00088       return exp(lnOfGamma(z)+lnOfGamma(w)-lnOfGamma(z+w));
00089     }
00090 
00091   template <class T> 
00092     T gammaP(T a, T x)
00093     {
00094       T gamser, gammcf, gln ;
00095       if(x < 0 || a <= 0){
00096 #ifdef USE_EXCEPTION       
00097         throw MatrixInputError();
00098 #else
00099         Error error("gammaP");
00100         error << "Invalid arguments to the incomplete gamma function.\n" ;
00101         error.fatal();
00102 #endif
00103       }
00104       if(x < (a+1)){
00105         return gammaSerie(a,x,gln);
00106       }
00107       return 1.0-gammaSerieCF(a,x,gln);
00108     }
00109 
00110   template <class T>
00111     T gammaQ(T a, T x)
00112     {
00113       T gamser, gammcf, gln ;
00114       if(x < 0 || a <= 0){
00115 #ifdef USE_EXCEPTION       
00116         throw MatrixInputError();
00117 #else
00118         Error error("gammaQ");
00119         error << "Invalid arguments to the incomplete gamma function.\n" ;
00120         error.fatal();
00121 #endif
00122       }
00123       if(x < (a+1)){
00124         return 1.0-gammaSerie(a,x,gln);
00125       }
00126       return gammaSerieCF(a,x,gln);
00127     }
00128 
00129   template <class T>
00130     T gammaSerie(T a, T x, T& gln)
00131     {
00132       int n;
00133       T sum, del, ap;
00134       gln=lnOfGamma(a);
00135       if(x<0){
00136 #ifdef USE_EXCEPTION       
00137         throw MatrixInputError();
00138 #else
00139         Error error("gammaSerie");
00140         error << "x can't be negative.\n" ;
00141         error.fatal();
00142         return 0; 
00143 #endif
00144       }
00145       
00146       ap = a ;
00147       del= sum = T(1)/a;
00148       for(n=MaximumIterations;n>0;--n){
00149         ++ap;
00150         del *= x/ap ;
00151         sum += del ;
00152         if(absolute(del)< absolute(sum)*Precision<T>())
00153           return sum*exp(-x+a*log(x)-gln) ;
00154       }
00155 #ifdef USING_EXCEPTION
00156       throw MatrixErr();
00157 #else
00158       Error error("gammaSerie");
00159       error << "a too large or MaximumIterations too small.\n" ;
00160       error.fatal();
00161       return 0;
00162 #endif
00163     }
00164 
00165   template <class T>
00166     T gammaSerieCF(T a, T x, T& gln)
00167     {
00168       int i;
00169       T an,b,c,d,del,h;
00170 
00171       gln=lnOfGamma(a);
00172 
00173       b= x + 1.0 - a;
00174       c= T(1)/MinimumFloat<T>();
00175       d= T(1)/b ;
00176       h= d ;
00177       for(i=1;i<=MaximumIterations;++i){
00178         an = -T(i)*T(i-a);
00179         b += 2 ;
00180         d = an*d + b ;
00181         if(absolute(d) < MinimumFloat<T>() )
00182           d = MinimumFloat<T>();
00183         c = b + an/c ;
00184         if(absolute(c) < MinimumFloat<T>() )
00185           c = MinimumFloat<T>();
00186         d=T(1)/d ;
00187         del=d*c;
00188         h *= del ;
00189         if(absolute(del-1.0) < Precision<T>()) break ;
00190         
00191       }
00192 
00193       if(i>MaximumIterations){
00194 #ifdef USING_EXCEPTION
00195         throw MatrixErr();
00196 #else
00197         Error error("gammaSerie");
00198         error << "a too large or MaximumIterations too small.\n" ;
00199         error.fatal();
00200         return 0;
00201 #endif
00202       }
00203       return exp(-x+a*log(x)-gln)*h ;
00204     }
00205 
00206   template <class T>
00207     T errorFcn(T x){
00208     return x < 0 ? -gammaP<T>(0.5,x*x) : gammaP<T>(0.5,x*x) ;
00209   }
00210 
00211   template <class T>
00212     T errorFcnC(T x){
00213     return x < 0 ? 1+gammaP<T>(0.5,x*x) : gammaQ<T>(0.5,x*x) ; 
00214   }
00215 
00216   template <class T>
00217     T errorFcnChebyshevC(T x)
00218     {
00219       T t,z,ans ;
00220       z = absolute(x);
00221       t = T(1)/(T(1)+0.5*z) ;
00222       ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*0.09678418+
00223                 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
00224                 t*(-0.82215223+t*0.17087277))))))));
00225       return x>0 ? ans : T(2)-ans ;
00226                                                  
00227     }
00228 
00229   template <class T>
00230     void kendallTau(const BasicArray<T>& data1, const BasicArray<T>& data2, T &tau, T &z, T& prob)
00231     {
00232       unsigned int n2=0,n1=0,k,j ;
00233       int is=0;
00234       T svar, aa, a2,a1 ;
00235       int n = data1.n() ;
00236 
00237       for(j=0;j<n-1;++j)
00238         for(k=j+1;k<n;++k){
00239           a1= data1[j]-data1[k] ;
00240           a2= data2[j]-data2[k] ;
00241           aa=a1*a2;
00242           if(aa){
00243             ++n1 ; //no tie
00244             ++n2 ;
00245             aa > 0 ? ++is : --is ;
00246           }
00247           else{
00248             if(a1) ++n1 ; // extra 'x' event
00249             if(a2) ++n2 ; // extra 'y' event
00250           }
00251         }
00252       tau = is/(T(sqrt((double)n1)*sqrt((double)n2))) ;
00253       svar = T(4*n+10)/T(9*n*(n-1)) ;
00254       z = tau/T(sqrt(svar)) ;
00255       prob = errorFcnChebyshevC<T>(absolute(z)/1.4142136) ;
00256     }
00257 
00258   void kendallTau(const BasicArray<int>& data1, const BasicArray<int>& data2, float &tau, float &z, float& prob)
00259     {
00260       unsigned int n2=0,n1=0,k,j ;
00261       int is=0;
00262       float svar, aa, a2,a1 ;
00263       int n = data1.n() ;
00264 
00265       
00266       for(j=0;j<n-1;++j)
00267         for(k=j+1;k<n;++k){
00268           a1= data1[j]-data1[k] ;
00269           a2= data2[j]-data2[k] ;
00270           aa=a1*a2;
00271           if(aa){
00272             ++n1 ; //no tie
00273             ++n2 ;
00274             aa > 0 ? ++is : --is ;
00275           }
00276           else{
00277             if(a1) ++n1 ; // extra 'x' event
00278             if(a2) ++n2 ; // extra 'y' event
00279           }
00280         }
00281       tau = is/(float(sqrt((double)n1)*sqrt((double)n2))) ;
00282       svar = float(4*n+10)/float(9*n*(n-1)) ;
00283       z = tau/float(sqrt(svar)) ;
00284       prob = errorFcnChebyshevC<float>(absolute(z)/1.4142136) ;
00285     }
00286 
00287   void kendallTau(const BasicArray<int>& data1, const BasicArray<int>& data2, double &tau, double &z, double& prob)
00288     {
00289       unsigned int n2=0,n1=0,k,j ;
00290       int is=0;
00291       double svar, aa, a2,a1 ;
00292       int n = data1.n() ;
00293 
00294       
00295       for(j=0;j<n-1;++j)
00296         for(k=j+1;k<n;++k){
00297           a1= data1[j]-data1[k] ;
00298           a2= data2[j]-data2[k] ;
00299           aa=a1*a2;
00300           if(aa){
00301             ++n1 ; //no tie
00302             ++n2 ;
00303             aa > 0 ? ++is : --is ;
00304           }
00305           else{
00306             if(a1) ++n1 ; // extra 'x' event
00307             if(a2) ++n2 ; // extra 'y' event
00308           }
00309         }
00310       tau = is/(double(sqrt((double)n1)*sqrt((double)n2))) ;
00311       svar = double(4*n+10)/double(9*n*(n-1)) ;
00312       z = tau/double(sqrt(svar)) ;
00313       prob = errorFcnChebyshevC<double>(absolute(z)/1.4142136) ;
00314     }
00315 
00316 }
00317 
00318 
00319 #endif