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

Generated on Tue Feb 2 09:33:32 2010 for IT++ by Doxygen 1.6.2