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