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