IT++ Logo

misc_stat.cpp

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

Generated on Sat Apr 19 10:41:15 2008 for IT++ by Doxygen 1.5.5