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
Generated on Tue Feb 2 09:33:32 2010 for IT++ by Doxygen 1.6.2