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 frobeniuos norm of a matrix for s = "fro" 00134 double norm(const mat &m, const std::string &s) 00135 { 00136 it_assert(s == "fro", "norm(): Unrecognised norm"); 00137 return std::sqrt(sum(diag(transpose(m) * m))); 00138 } 00139 00140 // Calculate the frobeniuos norm of a matrix for s = "fro" 00141 double norm(const cmat &m, const std::string &s) 00142 { 00143 it_assert(s == "fro", "norm(): Unrecognised norm"); 00144 return std::sqrt(sum(real(diag(hermitian_transpose(m) * m)))); 00145 } 00146 00147 00148 double variance(const cvec &v) 00149 { 00150 int len = v.size(); 00151 double sq_sum = 0.0; 00152 std::complex<double> sum = 0.0; 00153 const std::complex<double> *p = v._data(); 00154 00155 for (int i = 0; i < len; i++, p++) { 00156 sum += *p; 00157 sq_sum += std::norm(*p); 00158 } 00159 00160 return (double)(sq_sum - std::norm(sum) / len) / (len - 1); 00161 } 00162 00163 double moment(const vec &x, const int r) 00164 { 00165 double m = mean(x), mr = 0; 00166 int n = x.size(); 00167 double temp; 00168 00169 switch (r) { 00170 case 1: 00171 for (int j = 0; j < n; j++) 00172 mr += (x(j) - m); 00173 break; 00174 case 2: 00175 for (int j = 0; j < n; j++) 00176 mr += (x(j) - m) * (x(j) - m); 00177 break; 00178 case 3: 00179 for (int j = 0; j < n; j++) 00180 mr += (x(j) - m) * (x(j) - m) * (x(j) - m); 00181 break; 00182 case 4: 00183 for (int j = 0; j < n; j++) { 00184 temp = (x(j) - m) * (x(j) - m); 00185 temp *= temp; 00186 mr += temp; 00187 } 00188 break; 00189 default: 00190 for (int j = 0; j < n; j++) 00191 mr += std::pow(x(j) - m, double(r)); 00192 break; 00193 } 00194 00195 return mr / n; 00196 } 00197 00198 00199 double skewness(const vec &x) 00200 { 00201 int n = x.size(); 00202 00203 double k2 = variance(x) * n / (n - 1); // 2nd k-statistic 00204 double k3 = moment(x, 3) * n * n / (n - 1) / (n - 2); //3rd k-statistic 00205 00206 return k3 / std::pow(k2, 3.0 / 2.0); 00207 } 00208 00209 double kurtosisexcess(const vec &x) 00210 { 00211 int n = x.size(); 00212 double m2 = variance(x); 00213 double m4 = moment(x, 4); 00214 00215 double k2 = m2 * n / (n - 1); // 2nd k-statistic 00216 double k4 = (m4 * (n + 1) - 3 * (n - 1) * m2 * m2) * n * n / (n - 1) / (n - 2) / (n - 3); //4th k-statistic 00217 00218 return k4 / (k2*k2); 00219 } 00220 00221 } // namespace itpp
Generated on Wed Feb 9 2011 13:47:17 for IT++ by Doxygen 1.7.3