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
Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4