diff --git a/src/mlpack/core.hpp b/src/mlpack/core.hpp index d29a174724..a2275ad188 100644 --- a/src/mlpack/core.hpp +++ b/src/mlpack/core.hpp @@ -83,6 +83,7 @@ #include #include #include +#include #include #include #include diff --git a/src/mlpack/core/math/CMakeLists.txt b/src/mlpack/core/math/CMakeLists.txt index 6134d5f583..28131b2f31 100644 --- a/src/mlpack/core/math/CMakeLists.txt +++ b/src/mlpack/core/math/CMakeLists.txt @@ -13,6 +13,7 @@ set(SOURCES make_alias.hpp multiply_slices_impl.hpp multiply_slices.hpp + quantile.hpp random.hpp random.cpp random_basis.hpp diff --git a/src/mlpack/core/math/quantile.hpp b/src/mlpack/core/math/quantile.hpp new file mode 100644 index 0000000000..fc469257ab --- /dev/null +++ b/src/mlpack/core/math/quantile.hpp @@ -0,0 +1,123 @@ +/** + * @file core/math/quantile.hpp + * @author Shubham Agrawal + * + * Miscellaneous math quantile-related routines. + * + * mlpack is free software; you may redistribute it and/or modify it under the + * terms of the 3-clause BSD license. You should have received a copy of the + * 3-clause BSD license along with mlpack. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ +#ifndef MLPACK_CORE_MATH_QUANTILE_HPP +#define MLPACK_CORE_MATH_QUANTILE_HPP + +#include + +namespace mlpack { +namespace math /** Miscellaneous math routines. */ { + +/** + * Computes the inverse erf function using the rational approximation from + * Numerical Recipes. + * Code accompanying the article "Approximating the erfinv function" in + * GPU Computing Gems, Volume 2. + * + * @param x Input value. + */ +inline double ErfInverse(double x) +{ + double w, p; + + w = -log((1.0 - x) * (1.0 + x)); + + if (w < 6.250000) + { + w = w - 3.125000; + p = -3.6444120640178196996e-21; + p = -1.685059138182016589e-19 + p * w; + p = 1.2858480715256400167e-18 + p * w; + p = 1.115787767802518096e-17 + p * w; + p = -1.333171662854620906e-16 + p * w; + p = 2.0972767875968561637e-17 + p * w; + p = 6.6376381343583238325e-15 + p * w; + p = -4.0545662729752068639e-14 + p * w; + p = -8.1519341976054721522e-14 + p * w; + p = 2.6335093153082322977e-12 + p * w; + p = -1.2975133253453532498e-11 + p * w; + p = -5.4154120542946279317e-11 + p * w; + p = 1.051212273321532285e-09 + p * w; + p = -4.1126339803469836976e-09 + p * w; + p = -2.9070369957882005086e-08 + p * w; + p = 4.2347877827932403518e-07 + p * w; + p = -1.3654692000834678645e-06 + p * w; + p = -1.3882523362786468719e-05 + p * w; + p = 0.0001867342080340571352 + p * w; + p = -0.00074070253416626697512 + p * w; + p = -0.0060336708714301490533 + p * w; + p = 0.24015818242558961693 + p * w; + p = 1.6536545626831027356 + p * w; + } + else if (w < 16.000000) + { + w = sqrt(w) - 3.250000; + p = 2.2137376921775787049e-09; + p = 9.0756561938885390979e-08 + p * w; + p = -2.7517406297064545428e-07 + p * w; + p = 1.8239629214389227755e-08 + p * w; + p = 1.5027403968909827627e-06 + p * w; + p = -4.013867526981545969e-06 + p * w; + p = 2.9234449089955446044e-06 + p * w; + p = 1.2475304481671778723e-05 + p * w; + p = -4.7318229009055733981e-05 + p * w; + p = 6.8284851459573175448e-05 + p * w; + p = 2.4031110387097893999e-05 + p * w; + p = -0.0003550375203628474796 + p * w; + p = 0.00095328937973738049703 + p * w; + p = -0.0016882755560235047313 + p * w; + p = 0.0024914420961078508066 + p * w; + p = -0.0037512085075692412107 + p * w; + p = 0.005370914553590063617 + p * w; + p = 1.0052589676941592334 + p * w; + p = 3.0838856104922207635 + p * w; + } + else + { + w = sqrt(w) - 5.000000; + p = -2.7109920616438573243e-11; + p = -2.5556418169965252055e-10 + p * w; + p = 1.5076572693500548083e-09 + p * w; + p = -3.7894654401267369937e-09 + p * w; + p = 7.6157012080783393804e-09 + p * w; + p = -1.4960026627149240478e-08 + p * w; + p = 2.9147953450901080826e-08 + p * w; + p = -6.7711997758452339498e-08 + p * w; + p = 2.2900482228026654717e-07 + p * w; + p = -9.9298272942317002539e-07 + p * w; + p = 4.5260625972231537039e-06 + p * w; + p = -1.9681778105531670567e-05 + p * w; + p = 7.5995277030017761139e-05 + p * w; + p = -0.00021503011930044477347 + p * w; + p = -0.00013871931833623122026 + p * w; + p = 1.0103004648645343977 + p * w; + p = 4.8499064014085844221 + p * w; + } + return p * x; +} + +/** + * Computes the quantile function of Guassian distribution at given probability. + * + * @param p Probability value. + * @param mu Mean of the distribution. (Default 0) + * @param sigma Standard deviation of the distribution. (Default 1) + */ +inline double Quantile(double p, double mu = 0.0, double sigma = 1.0) +{ + return mu + sigma * std::sqrt(2.0) * ErfInverse(2 * p - 1); +} + +} // namespace math +} // namespace mlpack + +#endif diff --git a/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp b/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp index 1231850e85..878d621b1d 100644 --- a/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp +++ b/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp @@ -12,7 +12,7 @@ #include "cosine_tree.hpp" #include -#include +#include namespace mlpack { namespace tree { @@ -477,8 +477,7 @@ double CosineTree::MonteCarloError(CosineTree* node, // Fit a normal distribution using the calculated statistics, and calculate a // lower bound on the magnitudes for the passed 'delta' parameter. - boost::math::normal dist(mu, sigma); - double lowerBound = boost::math::quantile(dist, delta); + double lowerBound = math::Quantile(delta, mu, sigma); // Upper bound on the subspace reconstruction error. node->L2Error(node->FrobNormSquared() - lowerBound); diff --git a/src/mlpack/methods/kde/kde_rules_impl.hpp b/src/mlpack/methods/kde/kde_rules_impl.hpp index 38a719d4ee..062e740efa 100644 --- a/src/mlpack/methods/kde/kde_rules_impl.hpp +++ b/src/mlpack/methods/kde/kde_rules_impl.hpp @@ -17,7 +17,7 @@ #include "kde_rules.hpp" // Used for Monte Carlo estimation. -#include +#include namespace mlpack { namespace kde { @@ -191,9 +191,7 @@ Score(const size_t queryIndex, TreeType& referenceNode) // Monte Carlo probabilistic estimation. // Calculate z using accumulated alpha if possible. const double alpha = depthAlpha + accumMCAlpha(queryIndex); - const boost::math::normal normalDist; - const double z = - std::abs(boost::math::quantile(normalDist, alpha / 2)); + const double z = std::abs(math::Quantile(alpha / 2.0)); // Auxiliary variables. arma::vec sample; @@ -400,9 +398,7 @@ Score(TreeType& queryNode, TreeType& referenceNode) // Monte Carlo probabilistic estimation. // Calculate z using accumulated alpha if possible. const double alpha = depthAlpha + queryStat.AccumAlpha(); - const boost::math::normal normalDist; - const double z = - std::abs(boost::math::quantile(normalDist, alpha / 2)); + const double z = std::abs(math::Quantile(alpha / 2)); // Auxiliary variables. arma::vec sample; diff --git a/src/mlpack/methods/preprocess/preprocess_describe_main.cpp b/src/mlpack/methods/preprocess/preprocess_describe_main.cpp index 04faf5fe9c..2859190167 100644 --- a/src/mlpack/methods/preprocess/preprocess_describe_main.cpp +++ b/src/mlpack/methods/preprocess/preprocess_describe_main.cpp @@ -19,14 +19,12 @@ #include -#include -#include +#include using namespace mlpack; using namespace mlpack::data; using namespace mlpack::util; using namespace std; -using namespace boost; // Program Name. BINDING_USER_NAME("Descriptive Statistics"); @@ -187,25 +185,13 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers) // Load the data. arma::mat& data = params.Get("input"); - // Generate boost format recipe. - const string widthPrecision("%-" + to_string(width) + "." + - to_string(precision)); - const string widthOnly("%-" + to_string(width) + "."); - string stringFormat = ""; - string numberFormat = ""; - - // We are going to print 11 different categories. - for (size_t i = 0; i < 11; ++i) - { - stringFormat += widthOnly + "s"; - numberFormat += widthPrecision + "f"; - } - timers.Start("statistics"); // Print the headers. - Log::Info << boost::format(stringFormat) - % "dim" % "var" % "mean" % "std" % "median" % "min" % "max" - % "range" % "skew" % "kurt" % "SE" << endl; + Log::Info << setw(width) << "dim" << setw(width) << "var" << setw(width) << + "mean" << setw(width) << "std" << setw(width) << setw(width) << setw(width) << + "median" << setw(width) << "min" << setw(width) << "max" << setw(width) << + "range" << setw(width) << "skew" << setw(width) << "kurt" << setw(width) << + "SE" << endl; // Lambda function to print out the results. auto PrintStatResults = [&](size_t dim, bool rowMajor) @@ -223,19 +209,17 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers) const double fStd = arma::stddev(feature, population); // Print statistics of the given dimension. - Log::Info << boost::format(numberFormat) - % dim - % arma::var(feature, population) - % fMean - % fStd - % arma::median(feature) - % fMin - % fMax - % (fMax - fMin) // range - % Skewness(feature, fStd, fMean, population) - % Kurtosis(feature, fStd, fMean, population) - % StandardError(feature.n_elem, fStd) - << endl; + Log::Info << setprecision(precision) << setw(width) << dim << + setw(width) << arma::var(feature, population) << + setw(width) << fMean << + setw(width) << fStd << + setw(width) << arma::median(feature) << + setw(width) << fMin << + setw(width) << fMax << + setw(width) << (fMax - fMin) << + setw(width) << Skewness(feature, fStd, fMean, population) << + setw(width) << Kurtosis(feature, fStd, fMean, population) << + setw(width) << StandardError(feature.n_elem, fStd) << endl; }; // If the user specified dimension, describe statistics of the given diff --git a/src/mlpack/tests/main_tests/emst_test.cpp b/src/mlpack/tests/main_tests/emst_test.cpp index 3cf3a8802a..afc6269210 100644 --- a/src/mlpack/tests/main_tests/emst_test.cpp +++ b/src/mlpack/tests/main_tests/emst_test.cpp @@ -18,8 +18,6 @@ #include "../catch.hpp" -#include - using namespace mlpack; BINDING_TEST_FIXTURE(EMSTTestFixture); @@ -106,10 +104,10 @@ TEST_CASE_METHOD(EMSTTestFixture, "EMSTFirstTwoOutputRowsIntegerTest", for (size_t i = 0; i < params.Get("output").n_cols; ++i) { REQUIRE(params.Get("output")(0, i) == - Approx(boost::math::iround(params.Get("output")(0, i))). + Approx(std::round(params.Get("output")(0, i))). epsilon(1e-7)); REQUIRE(params.Get("output")(1, i) == - Approx(boost::math::iround(params.Get("output")(1, i))). + Approx(std::round(params.Get("output")(1, i))). epsilon(1e-7)); } }