IT++ Logo Newcom Logo

stat.cpp

Go to the documentation of this file.
00001 
00033 #include <itpp/base/stat.h>
00034 #include <itpp/base/svd.h>
00035 
00036 
00037 namespace itpp {
00038 
00039   double mean(const vec &v)
00040   {
00041     return sum(v)/v.length();
00042   }
00043 
00044   std::complex<double> mean(const cvec &v)
00045   {
00046     return sum(v)/double(v.size());
00047   }
00048 
00049   double mean(const svec &v)
00050   {
00051     return (double)sum(v)/v.length();
00052   }
00053 
00054   double mean(const ivec &v)
00055   {
00056     return (double)sum(v)/v.length();
00057   }
00058 
00059   double mean(const mat &m)
00060   {
00061     return sum(sum(m))/(m.rows()*m.cols());
00062   }
00063 
00064   std::complex<double> mean(const cmat &m)
00065   {
00066     return sum(sum(m))/static_cast<std::complex<double> >(m.rows()*m.cols());
00067   }
00068 
00069   double mean(const smat &m)
00070   {
00071     return static_cast<double>(sum(sum(m)))/(m.rows()*m.cols());
00072   }
00073 
00074   double mean(const imat &m)
00075   {
00076     return static_cast<double>(sum(sum(m)))/(m.rows()*m.cols());
00077   }
00078 
00079 
00080   double norm(const cvec &v)
00081   {
00082     double E = 0.0;
00083     for (int i = 0; i < v.length(); i++)
00084       E += std::norm(v[i]);
00085     
00086     return std::sqrt(E);
00087   }
00088 
00089   double norm(const cvec &v, int p)
00090   {
00091     double E = 0.0;
00092     for (int i = 0; i < v.size(); i++)
00093       E += std::pow(std::norm(v[i]), p / 2.0); // Yes, 2.0 is correct!
00094 
00095     return std::pow(E, 1.0 / p);
00096   }
00097 
00098   double norm(const cvec &v, const std::string &s) {
00099     return norm(v, 2);
00100   }
00101 
00102   /* 
00103    * Calculate the p-norm of a real matrix
00104    * p = 1: max(svd(m))
00105    * p = 2: max(sum(abs(X)))
00106    */
00107   double norm(const mat &m, int p)
00108   {
00109     it_assert((p == 1) || (p == 2), 
00110               "norm(): Can only calculate a matrix norm of order 1 or 2");
00111 
00112     if (p == 1)
00113       return max(sum(abs(m)));
00114     else
00115       return max(svd(m));
00116   }
00117 
00118   /*
00119    * Calculate the p-norm of a complex matrix
00120    * p = 1: max(svd(m))
00121    * p = 2: max(sum(abs(X)))
00122    */
00123   double norm(const cmat &m, int p)
00124   {
00125     it_assert((p == 1) || (p == 2), 
00126               "norm(): Can only calculate a matrix norm of order 1 or 2");
00127 
00128     if (p == 1)
00129       return max(sum(abs(m)));
00130     else
00131       return max(svd(m));
00132   }
00133 
00134   // Calculate the frobeniuos norm of a matrix for s = "fro"
00135   double norm(const mat &m, const std::string &s)
00136   {
00137     it_assert(s == "fro", "norm(): Unrecognised norm");
00138     return std::sqrt(sum(diag(transpose(m) * m)));
00139   }
00140 
00141   // Calculate the frobeniuos norm of a matrix for s = "fro"
00142   double norm(const cmat &m, const std::string &s)
00143   {
00144     it_assert(s == "fro", "norm(): Unrecognised norm");
00145     return std::sqrt(sum(real(diag(hermitian_transpose(m) * m))));
00146   }
00147 
00148 
00149   double variance(const cvec &v)
00150   {
00151     int len = v.size();
00152     double sq_sum=0.0;
00153     std::complex<double> sum=0.0;
00154     const std::complex<double> *p=v._data();
00155 
00156     for (int i=0; i<len; i++, p++) {
00157       sum += *p;
00158       sq_sum += std::norm(*p);
00159     }
00160 
00161     return (double)(sq_sum - std::norm(sum)/len) / (len-1);
00162   }
00163 
00164   double moment(const vec &x, const int r)
00165   {
00166     double m = mean(x), mr=0;
00167     int n = x.size();
00168     double temp;
00169 
00170       switch (r) {
00171       case 1:
00172         for (int j=0; j<n; j++)
00173           mr += (x(j)-m);
00174         break;
00175       case 2:
00176         for (int j=0; j<n; j++)
00177           mr += (x(j)-m) * (x(j)-m);
00178         break;
00179       case 3:
00180         for (int j=0; j<n; j++)
00181           mr += (x(j)-m) * (x(j)-m) * (x(j)-m);
00182         break;
00183       case 4:
00184         for (int j=0; j<n; j++) {
00185           temp = (x(j)-m) * (x(j)-m);
00186           temp *= temp;
00187           mr += temp;
00188         }
00189         break;
00190       default:
00191         for (int j=0; j<n; j++)
00192           mr += std::pow(x(j)-m, double(r));
00193         break;
00194       }
00195 
00196     return mr/n;
00197   }
00198 
00199 
00200   double skewness(const vec &x)
00201   {
00202     int n = x.size();
00203 
00204     double k2 = variance(x)*n/(n-1); // 2nd k-statistic
00205     double k3 = moment(x, 3)*n*n/(n-1)/(n-2); //3rd k-statistic
00206 
00207     return k3/std::pow(k2, 3.0/2.0);
00208   }
00209 
00210   double kurtosisexcess(const vec &x)
00211   {
00212     int n = x.size();
00213     double m2 = variance(x);
00214     double m4 = moment(x, 4);
00215 
00216     double k2 = m2*n/(n-1); // 2nd k-statistic
00217     double k4 = (m4*(n+1) - 3*(n-1)*m2*m2)*n*n/(n-1)/(n-2)/(n-3); //4th k-statistic
00218 
00219     return k4/(k2*k2);
00220   }
00221 
00222 } // namespace itpp
SourceForge Logo

Generated on Thu Apr 19 14:43:43 2007 for IT++ by Doxygen 1.5.1