IT++ Logo
misc_stat.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/algebra/svd.h>
00030 #include <itpp/stat/misc_stat.h>
00031 
00032 
00033 namespace itpp
00034 {
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 &)
00096 {
00097   return norm(v, 2);
00098 }
00099 
00100 /*
00101  * Calculate the p-norm of a real matrix
00102  * p = 1: max(svd(m))
00103  * p = 2: max(sum(abs(X)))
00104  */
00105 double norm(const mat &m, int p)
00106 {
00107   it_assert((p == 1) || (p == 2),
00108             "norm(): Can only calculate a matrix norm of order 1 or 2");
00109 
00110   if (p == 1)
00111     return max(sum(abs(m)));
00112   else
00113     return max(svd(m));
00114 }
00115 
00116 /*
00117  * Calculate the p-norm of a complex matrix
00118  * p = 1: max(svd(m))
00119  * p = 2: max(sum(abs(X)))
00120  */
00121 double norm(const cmat &m, int p)
00122 {
00123   it_assert((p == 1) || (p == 2),
00124             "norm(): Can only calculate a matrix norm of order 1 or 2");
00125 
00126   if (p == 1)
00127     return max(sum(abs(m)));
00128   else
00129     return max(svd(m));
00130 }
00131 
00132 // Calculate the Frobenius norm of matrix m for s = "fro"
00133 double norm(const mat &m, const std::string &s)
00134 {
00135   it_assert(s == "fro", "norm(): Unrecognised norm");
00136   double E = 0.0;
00137   for (int r = 0; r < m.rows(); ++r) {
00138     for (int c = 0; c < m.cols(); ++c) {
00139       E += m(r, c) * m(r, c);
00140     }
00141   }
00142   return std::sqrt(E);
00143 }
00144 
00145 // Calculate the Frobenius norm of matrix m for s = "fro"
00146 double norm(const cmat &m, const std::string &s)
00147 {
00148   it_assert(s == "fro", "norm(): Unrecognised norm");
00149   double E = 0.0;
00150   for (int r = 0; r < m.rows(); ++r) {
00151     for (int c = 0; c < m.cols(); ++c) {
00152       E += std::norm(m(r, c));
00153     }
00154   }
00155   return std::sqrt(E);
00156 }
00157 
00158 
00159 double variance(const cvec &v)
00160 {
00161   int len = v.size();
00162   double sq_sum = 0.0;
00163   std::complex<double> sum = 0.0;
00164   const std::complex<double> *p = v._data();
00165 
00166   for (int i = 0; i < len; i++, p++) {
00167     sum += *p;
00168     sq_sum += std::norm(*p);
00169   }
00170 
00171   return (double)(sq_sum - std::norm(sum) / len) / (len - 1);
00172 }
00173 
00174 double moment(const vec &x, const int r)
00175 {
00176   double m = mean(x), mr = 0;
00177   int n = x.size();
00178   double temp;
00179 
00180   switch (r) {
00181   case 1:
00182     for (int j = 0; j < n; j++)
00183       mr += (x(j) - m);
00184     break;
00185   case 2:
00186     for (int j = 0; j < n; j++)
00187       mr += (x(j) - m) * (x(j) - m);
00188     break;
00189   case 3:
00190     for (int j = 0; j < n; j++)
00191       mr += (x(j) - m) * (x(j) - m) * (x(j) - m);
00192     break;
00193   case 4:
00194     for (int j = 0; j < n; j++) {
00195       temp = (x(j) - m) * (x(j) - m);
00196       temp *= temp;
00197       mr += temp;
00198     }
00199     break;
00200   default:
00201     for (int j = 0; j < n; j++)
00202       mr += std::pow(x(j) - m, double(r));
00203     break;
00204   }
00205 
00206   return mr / n;
00207 }
00208 
00209 
00210 double skewness(const vec &x)
00211 {
00212   int n = x.size();
00213 
00214   double k2 = variance(x) * n / (n - 1); // 2nd k-statistic
00215   double k3 = moment(x, 3) * n * n / (n - 1) / (n - 2); //3rd k-statistic
00216 
00217   return k3 / std::pow(k2, 3.0 / 2.0);
00218 }
00219 
00220 double kurtosisexcess(const vec &x)
00221 {
00222   int n = x.size();
00223   double m2 = variance(x);
00224   double m4 = moment(x, 4);
00225 
00226   double k2 = m2 * n / (n - 1); // 2nd k-statistic
00227   double k4 = (m4 * (n + 1) - 3 * (n - 1) * m2 * m2) * n * n / (n - 1) / (n - 2) / (n - 3); //4th k-statistic
00228 
00229   return k4 / (k2*k2);
00230 }
00231 
00232 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:33 for IT++ by Doxygen 1.7.4