MLPACK
1.0.4
|
00001 00023 #ifndef __MLPACK_METHODS_MOG_PHI_HPP 00024 #define __MLPACK_METHODS_MOG_PHI_HPP 00025 00026 #include <mlpack/core.hpp> 00027 00028 namespace mlpack { 00029 namespace gmm { 00030 00046 inline double phi(const double x, const double mean, const double var) 00047 { 00048 return exp(-1.0 * ((x - mean) * (x - mean) / (2 * var))) 00049 / sqrt(2 * M_PI * var); 00050 } 00051 00068 inline double phi(const arma::vec& x, 00069 const arma::vec& mean, 00070 const arma::mat& cov) 00071 { 00072 arma::vec diff = mean - x; 00073 00074 // Parentheses required for Armadillo 3.0.0 bug. 00075 arma::vec exponent = -0.5 * (trans(diff) * inv(cov) * diff); 00076 00077 // TODO: What if det(cov) < 0? 00078 return pow(2 * M_PI, (double) x.n_elem / -2.0) * pow(det(cov), -0.5) * 00079 exp(exponent[0]); 00080 } 00081 00094 inline double phi(const arma::vec& x, 00095 const arma::vec& mean, 00096 const arma::mat& cov, 00097 const std::vector<arma::mat>& d_cov, 00098 arma::vec& g_mean, 00099 arma::vec& g_cov) 00100 { 00101 // We don't call out to another version of the function to avoid inverting the 00102 // covariance matrix more than once. 00103 arma::mat cinv = inv(cov); 00104 00105 arma::vec diff = mean - x; 00106 // Parentheses required for Armadillo 3.0.0 bug. 00107 arma::vec exponent = -0.5 * (trans(diff) * inv(cov) * diff); 00108 00109 long double f = pow(2 * M_PI, (double) x.n_elem / 2) * pow(det(cov), -0.5) 00110 * exp(exponent[0]); 00111 00112 // Calculate the g_mean values; this is a (1 x dim) vector. 00113 arma::vec invDiff = cinv * diff; 00114 g_mean = f * invDiff; 00115 00116 // Calculate the g_cov values; this is a (1 x (dim * (dim + 1) / 2)) vector. 00117 for (size_t i = 0; i < d_cov.size(); i++) 00118 { 00119 arma::mat inv_d = cinv * d_cov[i]; 00120 00121 g_cov[i] = f * dot(d_cov[i] * invDiff, invDiff) + 00122 accu(inv_d.diag()) / 2; 00123 } 00124 00125 return f; 00126 } 00127 00138 inline void phi(const arma::mat& x, 00139 const arma::vec& mean, 00140 const arma::mat& cov, 00141 arma::vec& probabilities) 00142 { 00143 // Column i of 'diffs' is the difference between x.col(i) and the mean. 00144 arma::mat diffs = x - (mean * arma::ones<arma::rowvec>(x.n_cols)); 00145 00146 // Now, we only want to calculate the diagonal elements of (diffs' * cov^-1 * 00147 // diffs). We just don't need any of the other elements. We can calculate 00148 // the right hand part of the equation (instead of the left side) so that 00149 // later we are referencing columns, not rows -- that is faster. 00150 arma::mat rhs = -0.5 * inv(cov) * diffs; 00151 arma::vec exponents(diffs.n_cols); // We will now fill this. 00152 for (size_t i = 0; i < diffs.n_cols; i++) 00153 exponents(i) = exp(accu(diffs.unsafe_col(i) % rhs.unsafe_col(i))); 00154 00155 probabilities = pow(2 * M_PI, (double) mean.n_elem / -2.0) * 00156 pow(det(cov), -0.5) * exponents; 00157 } 00158 00159 }; // namespace gmm 00160 }; // namespace mlpack 00161 00162 #endif