diff --git a/include/core/CRapidJsonWriterBase.h b/include/core/CRapidJsonWriterBase.h
index 2e6ebc0019..71a3da09cd 100644
--- a/include/core/CRapidJsonWriterBase.h
+++ b/include/core/CRapidJsonWriterBase.h
@@ -20,10 +20,10 @@
 #include <rapidjson/writer.h>
 
 #include <boost/iterator/indirect_iterator.hpp>
-#include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/unordered/unordered_map.hpp>
 #include <boost/unordered/unordered_set.hpp>
 
+#include <cmath>
 #include <memory>
 #include <stack>
 
@@ -166,7 +166,7 @@ class CRapidJsonWriterBase
 
     bool Double(double d) {
         // rewrite NaN and Infinity to 0
-        if (!(boost::math::isfinite)(d)) {
+        if (std::isfinite(d) == false) {
             return TRapidJsonWriterBase::Int(0);
         }
 
@@ -355,7 +355,7 @@ class CRapidJsonWriterBase
     //! Adds a double field with the name fieldname to an object.
     //! \p fieldName must outlive \p obj or memory corruption will occur.
     void addDoubleFieldToObj(const std::string& fieldName, double value, TValue& obj) const {
-        if (!(boost::math::isfinite)(value)) {
+        if (std::isfinite(value) == false) {
             LOG_ERROR(<< "Adding " << value << " to the \"" << fieldName
                       << "\" field of a JSON document");
             // Don't return - make a best effort to add the value
@@ -489,7 +489,7 @@ class CRapidJsonWriterBase
     //! Log a message if we're trying to add nan/infinity to a JSON array
     template<typename NUMBER>
     void checkArrayNumberFinite(NUMBER val, const std::string& fieldName, bool& considerLogging) const {
-        if (considerLogging && !(boost::math::isfinite)(val)) {
+        if (considerLogging && (std::isfinite(val) == false)) {
             LOG_ERROR(<< "Adding " << val << " to the \"" << fieldName
                       << "\" array in a JSON document");
             // Don't return - make a best effort to add the value
diff --git a/include/maths/CMathsFuncs.h b/include/maths/CMathsFuncs.h
index 6b5833fc76..f0ae0aa434 100644
--- a/include/maths/CMathsFuncs.h
+++ b/include/maths/CMathsFuncs.h
@@ -38,8 +38,7 @@ namespace maths {
 //!
 class MATHS_EXPORT CMathsFuncs : private core::CNonInstantiatable {
 public:
-    //! Wrapper around boost::math::isnan() which avoids the need to add
-    //! cryptic brackets everywhere to deal with macros.
+    //! Check if value is NaN.
     static bool isNan(double val);
     //! Check if any of the components are NaN.
     template<std::size_t N>
@@ -54,8 +53,7 @@ class MATHS_EXPORT CMathsFuncs : private core::CNonInstantiatable {
     //! Check if an element is NaN.
     static bool isNan(const core::CSmallVectorBase<double>& val);
 
-    //! Wrapper around boost::math::isinf() which avoids the need to add
-    //! cryptic brackets everywhere to deal with macros.
+    //! Check if value is infinite.
     static bool isInf(double val);
     //! Check if any of the components are infinite.
     template<std::size_t N>
diff --git a/include/maths/CMultivariateNormalConjugate.h b/include/maths/CMultivariateNormalConjugate.h
index 6c28b77d48..07007717d5 100644
--- a/include/maths/CMultivariateNormalConjugate.h
+++ b/include/maths/CMultivariateNormalConjugate.h
@@ -28,11 +28,11 @@
 #include <maths/CNormalMeanPrecConjugate.h>
 #include <maths/CRestoreParams.h>
 #include <maths/CSampling.h>
+#include <maths/CTools.h>
 #include <maths/ProbabilityAggregators.h>
 
 #include <boost/make_unique.hpp>
 #include <boost/math/special_functions/beta.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 #include <boost/optional.hpp>
 
@@ -1079,41 +1079,42 @@ class CMultivariateNormalConjugate : public CMultivariatePrior {
             LOG_ERROR(<< "Failed to calculate log det " << wishartScaleMatrixPost);
             return maths_t::E_FpFailed;
         }
+        double logGammaPost = 0.0;
+        double logGammaPrior = 0.0;
+        double logGammaPostMinusPrior = 0.0;
 
-        try {
-            double logGammaPostMinusPrior = 0.0;
-            for (std::size_t i = 0u; i < N; ++i) {
-                logGammaPostMinusPrior +=
-                    boost::math::lgamma(
-                        0.5 * (wishartDegreesFreedomPost - static_cast<double>(i))) -
-                    boost::math::lgamma(0.5 * (wishartDegreesFreedomPrior -
-                                               static_cast<double>(i)));
+        for (std::size_t i = 0u; i < N; ++i) {
+            if (CTools::lgamma(0.5 * (wishartDegreesFreedomPost - static_cast<double>(i)),
+                               logGammaPost, true) &&
+                CTools::lgamma(0.5 * (wishartDegreesFreedomPrior - static_cast<double>(i)),
+                               logGammaPrior, true)) {
+
+                logGammaPostMinusPrior += logGammaPost - logGammaPrior;
+            } else {
+                return maths_t::E_FpOverflowed;
             }
-            LOG_TRACE(<< "numberSamples = " << numberSamples);
-            LOG_TRACE(<< "logGaussianPrecisionPrior = " << logGaussianPrecisionPrior
-                      << ", logGaussianPrecisionPost  = " << logGaussianPrecisionPost);
-            LOG_TRACE(<< "wishartDegreesFreedomPrior = " << wishartDegreesFreedomPrior
-                      << ", wishartDegreesFreedomPost  = " << wishartDegreesFreedomPost);
-            LOG_TRACE(<< "wishartScaleMatrixPrior = " << m_WishartScaleMatrix);
-            LOG_TRACE(<< "wishartScaleMatrixPost  = " << wishartScaleMatrixPost);
-            LOG_TRACE(<< "logDeterminantPrior = " << logDeterminantPrior
-                      << ", logDeterminantPost = " << logDeterminantPost);
-            LOG_TRACE(<< "logGammaPostMinusPrior = " << logGammaPostMinusPrior);
-            LOG_TRACE(<< "logCountVarianceScales = " << logCountVarianceScales);
-
-            double d = static_cast<double>(N);
-            result = 0.5 * (wishartDegreesFreedomPrior * logDeterminantPrior -
-                            wishartDegreesFreedomPost * logDeterminantPost -
-                            d * (logGaussianPrecisionPost - logGaussianPrecisionPrior) +
-                            (wishartDegreesFreedomPost - wishartDegreesFreedomPrior) *
-                                d * core::constants::LOG_TWO +
-                            2.0 * logGammaPostMinusPrior -
-                            numberSamples * d * core::constants::LOG_TWO_PI -
-                            logCountVarianceScales);
-        } catch (const std::exception& e) {
-            LOG_ERROR(<< "Failed to calculate marginal likelihood: " << e.what());
-            return maths_t::E_FpFailed;
         }
+        LOG_TRACE(<< "numberSamples = " << numberSamples);
+        LOG_TRACE(<< "logGaussianPrecisionPrior = " << logGaussianPrecisionPrior
+                  << ", logGaussianPrecisionPost  = " << logGaussianPrecisionPost);
+        LOG_TRACE(<< "wishartDegreesFreedomPrior = " << wishartDegreesFreedomPrior
+                  << ", wishartDegreesFreedomPost  = " << wishartDegreesFreedomPost);
+        LOG_TRACE(<< "wishartScaleMatrixPrior = " << m_WishartScaleMatrix);
+        LOG_TRACE(<< "wishartScaleMatrixPost  = " << wishartScaleMatrixPost);
+        LOG_TRACE(<< "logDeterminantPrior = " << logDeterminantPrior
+                  << ", logDeterminantPost = " << logDeterminantPost);
+        LOG_TRACE(<< "logGammaPostMinusPrior = " << logGammaPostMinusPrior);
+        LOG_TRACE(<< "logCountVarianceScales = " << logCountVarianceScales);
+
+        double d = static_cast<double>(N);
+        result = 0.5 * (wishartDegreesFreedomPrior * logDeterminantPrior -
+                        wishartDegreesFreedomPost * logDeterminantPost -
+                        d * (logGaussianPrecisionPost - logGaussianPrecisionPrior) +
+                        (wishartDegreesFreedomPost - wishartDegreesFreedomPrior) *
+                            d * core::constants::LOG_TWO +
+                        2.0 * logGammaPostMinusPrior -
+                        numberSamples * d * core::constants::LOG_TWO_PI - logCountVarianceScales);
+
         return static_cast<maths_t::EFloatingPointErrorStatus>(CMathsFuncs::fpStatus(result));
     }
 
diff --git a/include/maths/CTools.h b/include/maths/CTools.h
index 890120b566..60e726e8cf 100644
--- a/include/maths/CTools.h
+++ b/include/maths/CTools.h
@@ -695,7 +695,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
     static double logOneMinusX(double x);
 
     //! A wrapper around lgamma which handles corner cases if requested
-    static bool lgamma(double value, double& result, bool checkForFinite = false);
+    static bool lgamma(double value, double& result, bool checkForFinite = true);
 };
 }
 }
diff --git a/lib/maths/CCategoricalTools.cc b/lib/maths/CCategoricalTools.cc
index 6f38196b73..9a238893a5 100644
--- a/lib/maths/CCategoricalTools.cc
+++ b/lib/maths/CCategoricalTools.cc
@@ -14,7 +14,6 @@
 #include <maths/MathsTypes.h>
 
 #include <boost/math/distributions/binomial.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 
 #include <cmath>
@@ -408,8 +407,7 @@ double CCategoricalTools::logBinomialCoefficient(std::size_t n, std::size_t m) {
     }
     double n_ = static_cast<double>(n);
     double m_ = static_cast<double>(m);
-    return boost::math::lgamma(n_ + 1.0) - boost::math::lgamma(m_ + 1.0) -
-           boost::math::lgamma(n_ - m_ + 1.0);
+    return std::lgamma(n_ + 1.0) - std::lgamma(m_ + 1.0) - std::lgamma(n_ - m_ + 1.0);
 }
 
 double CCategoricalTools::binomialCoefficient(std::size_t n, std::size_t m) {
@@ -544,8 +542,19 @@ CCategoricalTools::logBinomialProbability(std::size_t n, double p, std::size_t m
 
     double n_ = static_cast<double>(n);
     double m_ = static_cast<double>(m);
-    result = std::min(boost::math::lgamma(n_ + 1.0) - boost::math::lgamma(m_ + 1.0) -
-                          boost::math::lgamma(n_ - m_ + 1.0) +
+
+    double logGammaNPlusOne = 0.0;
+    double logGammaMPlusOne = 0.0;
+    double logGammaNMinusMPlusOne = 0.0;
+
+    if ((CTools::lgamma(n_ + 1.0, logGammaNPlusOne, true) &&
+         CTools::lgamma(m_ + 1.0, logGammaMPlusOne, true) &&
+         CTools::lgamma(n_ - m_ + 1.0, logGammaNMinusMPlusOne, true)) == false) {
+
+        return maths_t::E_FpOverflowed;
+    }
+
+    result = std::min(logGammaNPlusOne - logGammaMPlusOne - logGammaNMinusMPlusOne +
                           m_ * std::log(p) + (n_ - m_) * std::log(1.0 - p),
                       0.0);
     return maths_t::E_FpNoErrors;
@@ -570,7 +579,11 @@ CCategoricalTools::logMultinomialProbability(const TDoubleVec& probabilities,
     }
 
     double n_ = static_cast<double>(n);
-    double logP = boost::math::lgamma(n_ + 1.0);
+    double logP = 0.0;
+
+    if (CTools::lgamma(n_ + 1.0, logP, true) == false) {
+        return maths_t::E_FpOverflowed;
+    }
 
     for (std::size_t i = 0u; i < ni.size(); ++i) {
         double ni_ = static_cast<double>(ni[i]);
@@ -584,7 +597,13 @@ CCategoricalTools::logMultinomialProbability(const TDoubleVec& probabilities,
                 result = boost::numeric::bounds<double>::lowest();
                 return maths_t::E_FpOverflowed;
             }
-            logP += ni_ * std::log(pi_) - boost::math::lgamma(ni_ + 1.0);
+
+            double logGammaNiPlusOne = 0.0;
+            if (CTools::lgamma(ni_ + 1.0, logGammaNiPlusOne, true) == false) {
+                return maths_t::E_FpOverflowed;
+            }
+
+            logP += ni_ * std::log(pi_) - logGammaNiPlusOne;
         }
     }
 
diff --git a/lib/maths/CGammaRateConjugate.cc b/lib/maths/CGammaRateConjugate.cc
index 2e49c1a825..78cf2bfc06 100644
--- a/lib/maths/CGammaRateConjugate.cc
+++ b/lib/maths/CGammaRateConjugate.cc
@@ -27,7 +27,6 @@
 #include <boost/math/distributions/beta.hpp>
 #include <boost/math/distributions/gamma.hpp>
 #include <boost/math/special_functions/digamma.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/tools/roots.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 
@@ -637,35 +636,37 @@ class CLogMarginalLikelihood : core::CNonCopyable {
         double logGammaScaledLikelihoodShape = 0.0;
         double scaledImpliedShape = 0.0;
 
-        try {
-            for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
-                double n = maths_t::countForUpdate(m_Weights[i]);
-                double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
-                                       maths_t::countVarianceScale(m_Weights[i]);
-                m_NumberSamples += n;
-                if (varianceScale != 1.0) {
-                    logVarianceScaleSum -= m_LikelihoodShape / varianceScale *
-                                           std::log(varianceScale);
-                    logGammaScaledLikelihoodShape +=
-                        n * boost::math::lgamma(m_LikelihoodShape / varianceScale);
-                    scaledImpliedShape += n * m_LikelihoodShape / varianceScale;
-                } else {
-                    nResidual += n;
-                }
+        for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
+            double n = maths_t::countForUpdate(m_Weights[i]);
+            double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
+                                   maths_t::countVarianceScale(m_Weights[i]);
+            m_NumberSamples += n;
+            if (varianceScale != 1.0) {
+                logVarianceScaleSum -= m_LikelihoodShape / varianceScale *
+                                       std::log(varianceScale);
+                logGammaScaledLikelihoodShape +=
+                    n * std::lgamma(m_LikelihoodShape / varianceScale);
+                scaledImpliedShape += n * m_LikelihoodShape / varianceScale;
+            } else {
+                nResidual += n;
             }
+        }
 
-            m_ImpliedShape = scaledImpliedShape + nResidual * m_LikelihoodShape + m_PriorShape;
+        m_ImpliedShape = scaledImpliedShape + nResidual * m_LikelihoodShape + m_PriorShape;
 
-            LOG_TRACE(<< "numberSamples = " << m_NumberSamples);
+        LOG_TRACE(<< "numberSamples = " << m_NumberSamples);
 
-            m_Constant = m_PriorShape * std::log(m_PriorRate) -
-                         boost::math::lgamma(m_PriorShape) +
-                         logVarianceScaleSum - logGammaScaledLikelihoodShape -
-                         nResidual * boost::math::lgamma(m_LikelihoodShape) +
-                         boost::math::lgamma(m_ImpliedShape);
-        } catch (const std::exception& e) {
-            LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
+        m_Constant = m_PriorShape * std::log(m_PriorRate) - std::lgamma(m_PriorShape) +
+                     logVarianceScaleSum - logGammaScaledLikelihoodShape -
+                     nResidual * std::lgamma(m_LikelihoodShape) +
+                     std::lgamma(m_ImpliedShape);
+
+        if (std::isnan(m_ImpliedShape) || std::isnan(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood: floating point nan");
             this->addErrorStatus(maths_t::E_FpFailed);
+        } else if (std::isinf(m_ImpliedShape) || std::isinf(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood: floating point overflow");
+            this->addErrorStatus(maths_t::E_FpOverflowed);
         }
     }
 
diff --git a/lib/maths/CLogNormalMeanPrecConjugate.cc b/lib/maths/CLogNormalMeanPrecConjugate.cc
index a4d8271788..3864113fd0 100644
--- a/lib/maths/CLogNormalMeanPrecConjugate.cc
+++ b/lib/maths/CLogNormalMeanPrecConjugate.cc
@@ -30,7 +30,6 @@
 #include <boost/math/distributions/lognormal.hpp>
 #include <boost/math/distributions/normal.hpp>
 #include <boost/math/distributions/students_t.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 
 #include <algorithm>
@@ -487,49 +486,51 @@ class CLogMarginalLikelihood : core::CNonCopyable {
 private:
     //! Compute all the constants in the integrand.
     void precompute() {
-        try {
-            double logVarianceScaleSum = 0.0;
-
-            if (maths_t::hasSeasonalVarianceScale(m_Weights) ||
-                maths_t::hasCountVarianceScale(m_Weights)) {
-                m_Scales.reserve(m_Weights.size());
-                double r = m_Rate / m_Shape;
-                double s = std::exp(-r);
-                for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
-                    double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
-                                           maths_t::countVarianceScale(m_Weights[i]);
-
-                    // Get the scale and shift of the exponentiated Gaussian.
-                    if (varianceScale == 1.0) {
-                        m_Scales.emplace_back(1.0, 0.0);
-                    } else {
-                        double t = r + std::log(s + varianceScale * (1.0 - s));
-                        m_Scales.emplace_back(t / r, 0.5 * (r - t));
-                        logVarianceScaleSum += std::log(t / r);
-                    }
+        double logVarianceScaleSum = 0.0;
+
+        if (maths_t::hasSeasonalVarianceScale(m_Weights) ||
+            maths_t::hasCountVarianceScale(m_Weights)) {
+            m_Scales.reserve(m_Weights.size());
+            double r = m_Rate / m_Shape;
+            double s = std::exp(-r);
+            for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
+                double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
+                                       maths_t::countVarianceScale(m_Weights[i]);
+
+                // Get the scale and shift of the exponentiated Gaussian.
+                if (varianceScale == 1.0) {
+                    m_Scales.emplace_back(1.0, 0.0);
+                } else {
+                    double t = r + std::log(s + varianceScale * (1.0 - s));
+                    m_Scales.emplace_back(t / r, 0.5 * (r - t));
+                    logVarianceScaleSum += std::log(t / r);
                 }
             }
+        }
 
-            m_NumberSamples = 0.0;
-            double weightedNumberSamples = 0.0;
+        m_NumberSamples = 0.0;
+        double weightedNumberSamples = 0.0;
 
-            for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
-                double n = maths_t::countForUpdate(m_Weights[i]);
-                m_NumberSamples += n;
-                weightedNumberSamples +=
-                    n / (m_Scales.empty() ? 1.0 : m_Scales[i].first);
-            }
+        for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
+            double n = maths_t::countForUpdate(m_Weights[i]);
+            m_NumberSamples += n;
+            weightedNumberSamples += n / (m_Scales.empty() ? 1.0 : m_Scales[i].first);
+        }
 
-            double impliedShape = m_Shape + 0.5 * m_NumberSamples;
-            double impliedPrecision = m_Precision + weightedNumberSamples;
+        double impliedShape = m_Shape + 0.5 * m_NumberSamples;
+        double impliedPrecision = m_Precision + weightedNumberSamples;
 
-            m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
-                         0.5 * m_NumberSamples * LOG_2_PI - 0.5 * logVarianceScaleSum +
-                         boost::math::lgamma(impliedShape) -
-                         boost::math::lgamma(m_Shape) + m_Shape * std::log(m_Rate);
-        } catch (const std::exception& e) {
-            LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
+        m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
+                     0.5 * m_NumberSamples * LOG_2_PI -
+                     0.5 * logVarianceScaleSum + std::lgamma(impliedShape) -
+                     std::lgamma(m_Shape) + m_Shape * std::log(m_Rate);
+
+        if (std::isnan(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood, floating point nan");
             this->addErrorStatus(maths_t::E_FpFailed);
+        } else if (std::isinf(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood, floating point overflow");
+            this->addErrorStatus(maths_t::E_FpOverflowed);
         }
     }
 
@@ -1219,7 +1220,7 @@ void CLogNormalMeanPrecConjugate::sampleMarginalLikelihood(std::size_t numberSam
             double z = (xq - m_GaussianMean - scale * scale) / scale /
                        boost::math::double_constants::root_two;
 
-            double partialExpectation = mean * (1.0 + boost::math::erf(z)) / 2.0;
+            double partialExpectation = mean * (1.0 + std::erf(z)) / 2.0;
 
             double sample = static_cast<double>(numberSamples) *
                                 (partialExpectation - lastPartialExpectation) -
diff --git a/lib/maths/CMathsFuncs.cc b/lib/maths/CMathsFuncs.cc
index 9df66f3d48..1a91614c08 100644
--- a/lib/maths/CMathsFuncs.cc
+++ b/lib/maths/CMathsFuncs.cc
@@ -8,21 +8,13 @@
 
 #include <maths/CMathsFuncsForMatrixAndVectorTypes.h>
 
-#include <boost/math/special_functions/fpclassify.hpp>
-
-#ifdef isnan
-#undef isnan
-#endif
-
-#ifdef isinf
-#undef isinf
-#endif
+#include <cmath>
 
 namespace ml {
 namespace maths {
 
 bool CMathsFuncs::isNan(double val) {
-    return boost::math::isnan(val);
+    return std::isnan(val);
 }
 bool CMathsFuncs::isNan(const CSymmetricMatrix<double>& val) {
     return anElement(static_cast<bool (*)(double)>(&isNan), val);
@@ -40,7 +32,7 @@ bool CMathsFuncs::isNan(const core::CSmallVectorBase<double>& val) {
 }
 
 bool CMathsFuncs::isInf(double val) {
-    return boost::math::isinf(val);
+    return std::isinf(val);
 }
 bool CMathsFuncs::isInf(const CVector<double>& val) {
     return aComponent(static_cast<bool (*)(double)>(&isInf), val);
@@ -58,7 +50,7 @@ bool CMathsFuncs::isInf(const core::CSmallVectorBase<double>& val) {
 }
 
 bool CMathsFuncs::isFinite(double val) {
-    return boost::math::isfinite(val);
+    return std::isfinite(val);
 }
 bool CMathsFuncs::isFinite(const CVector<double>& val) {
     return everyComponent(static_cast<bool (*)(double)>(&isFinite), val);
diff --git a/lib/maths/CMultinomialConjugate.cc b/lib/maths/CMultinomialConjugate.cc
index c220409e17..01bf4a2d5e 100644
--- a/lib/maths/CMultinomialConjugate.cc
+++ b/lib/maths/CMultinomialConjugate.cc
@@ -25,7 +25,6 @@
 #include <maths/ProbabilityAggregators.h>
 
 #include <boost/math/distributions/beta.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 #include <boost/range.hpp>
 #include <boost/tuple/tuple.hpp>
@@ -656,37 +655,28 @@ CMultinomialConjugate::jointLogMarginalLikelihood(const TDouble1Vec& samples,
         categoryCounts[samples[i]] += n;
     }
 
-    try {
-        LOG_TRACE(<< "# samples = " << numberSamples
-                  << ", total concentration = " << m_TotalConcentration);
+    LOG_TRACE(<< "# samples = " << numberSamples
+              << ", total concentration = " << m_TotalConcentration);
 
-        result = boost::math::lgamma(numberSamples + 1.0) +
-                 boost::math::lgamma(m_TotalConcentration) -
-                 boost::math::lgamma(m_TotalConcentration + numberSamples);
+    result = std::lgamma(numberSamples + 1.0) + std::lgamma(m_TotalConcentration) -
+             std::lgamma(m_TotalConcentration + numberSamples);
 
-        for (TDoubleDoubleMapCItr countItr = categoryCounts.begin();
-             countItr != categoryCounts.end(); ++countItr) {
-            double category = countItr->first;
-            double count = countItr->second;
-            LOG_TRACE(<< "category = " << category << ", count = " << count);
+    for (TDoubleDoubleMapCItr countItr = categoryCounts.begin();
+         countItr != categoryCounts.end(); ++countItr) {
+        double category = countItr->first;
+        double count = countItr->second;
+        LOG_TRACE(<< "category = " << category << ", count = " << count);
 
-            result -= boost::math::lgamma(countItr->second + 1.0);
+        result -= std::lgamma(countItr->second + 1.0);
 
-            std::size_t index = std::lower_bound(m_Categories.begin(),
-                                                 m_Categories.end(), category) -
-                                m_Categories.begin();
-            if (index < m_Categories.size() && m_Categories[index] == category) {
-                LOG_TRACE(<< "concentration = " << m_Concentrations[index]);
-                result += boost::math::lgamma(m_Concentrations[index] + count) -
-                          boost::math::lgamma(m_Concentrations[index]);
-            }
+        std::size_t index = std::lower_bound(m_Categories.begin(),
+                                             m_Categories.end(), category) -
+                            m_Categories.begin();
+        if (index < m_Categories.size() && m_Categories[index] == category) {
+            LOG_TRACE(<< "concentration = " << m_Concentrations[index]);
+            result += std::lgamma(m_Concentrations[index] + count) -
+                      std::lgamma(m_Concentrations[index]);
         }
-    } catch (const std::exception& e) {
-        LOG_ERROR(<< "Unable to compute joint log likelihood: " << e.what()
-                  << ", samples = " << core::CContainerPrinter::print(samples)
-                  << ", categories = " << core::CContainerPrinter::print(m_Categories)
-                  << ", concentrations = " << core::CContainerPrinter::print(m_Concentrations));
-        return maths_t::E_FpFailed;
     }
 
     LOG_TRACE(<< "result = " << result);
diff --git a/lib/maths/CNormalMeanPrecConjugate.cc b/lib/maths/CNormalMeanPrecConjugate.cc
index 1dd2332e45..b81b7fb4f1 100644
--- a/lib/maths/CNormalMeanPrecConjugate.cc
+++ b/lib/maths/CNormalMeanPrecConjugate.cc
@@ -356,40 +356,41 @@ class CLogMarginalLikelihood : core::CNonCopyable {
         TMeanVarAccumulator sampleMoments;
         double logVarianceScaleSum = 0.0;
 
-        try {
-            for (std::size_t i = 0u; i < samples.size(); ++i) {
-                double n = maths_t::countForUpdate(weights[i]);
-                double seasonalScale =
-                    std::sqrt(maths_t::seasonalVarianceScale(weights[i]));
-                double countVarianceScale = maths_t::countVarianceScale(weights[i]);
-                double w = 1.0 / countVarianceScale;
-                m_NumberSamples += n;
-                if (seasonalScale != 1.0) {
-                    sampleMoments.add(predictionMean + (samples[i] - predictionMean) / seasonalScale,
-                                      n * w);
-                    logVarianceScaleSum += 2.0 * std::log(seasonalScale);
-                } else {
-                    sampleMoments.add(samples[i], n * w);
-                }
-                if (countVarianceScale != 1.0) {
-                    logVarianceScaleSum += std::log(countVarianceScale);
-                }
+        for (std::size_t i = 0u; i < samples.size(); ++i) {
+            double n = maths_t::countForUpdate(weights[i]);
+            double seasonalScale = std::sqrt(maths_t::seasonalVarianceScale(weights[i]));
+            double countVarianceScale = maths_t::countVarianceScale(weights[i]);
+            double w = 1.0 / countVarianceScale;
+            m_NumberSamples += n;
+            if (seasonalScale != 1.0) {
+                sampleMoments.add(predictionMean + (samples[i] - predictionMean) / seasonalScale,
+                                  n * w);
+                logVarianceScaleSum += 2.0 * std::log(seasonalScale);
+            } else {
+                sampleMoments.add(samples[i], n * w);
             }
-            m_WeightedNumberSamples = CBasicStatistics::count(sampleMoments);
-            m_SampleMean = CBasicStatistics::mean(sampleMoments);
-            m_SampleSquareDeviation = (m_WeightedNumberSamples - 1.0) *
-                                      CBasicStatistics::variance(sampleMoments);
-
-            double impliedShape = m_Shape + 0.5 * m_NumberSamples;
-            double impliedPrecision = m_Precision + m_WeightedNumberSamples;
-
-            m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
-                         0.5 * m_NumberSamples * LOG_2_PI - 0.5 * logVarianceScaleSum +
-                         boost::math::lgamma(impliedShape) -
-                         boost::math::lgamma(m_Shape) + m_Shape * std::log(m_Rate);
-        } catch (const std::exception& e) {
-            LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
+            if (countVarianceScale != 1.0) {
+                logVarianceScaleSum += std::log(countVarianceScale);
+            }
+        }
+        m_WeightedNumberSamples = CBasicStatistics::count(sampleMoments);
+        m_SampleMean = CBasicStatistics::mean(sampleMoments);
+        m_SampleSquareDeviation = (m_WeightedNumberSamples - 1.0) *
+                                  CBasicStatistics::variance(sampleMoments);
+
+        double impliedShape = m_Shape + 0.5 * m_NumberSamples;
+        double impliedPrecision = m_Precision + m_WeightedNumberSamples;
+
+        m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
+                     0.5 * m_NumberSamples * LOG_2_PI -
+                     0.5 * logVarianceScaleSum + std::lgamma(impliedShape) -
+                     std::lgamma(m_Shape) + m_Shape * std::log(m_Rate);
+        if (std::isnan(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood, floating point nan");
             this->addErrorStatus(maths_t::E_FpFailed);
+        } else if (std::isinf(m_Constant)) {
+            LOG_ERROR(<< "Error calculating marginal likelihood, floating point overflow");
+            this->addErrorStatus(maths_t::E_FpOverflowed);
         }
     }
 
diff --git a/lib/maths/CPoissonMeanConjugate.cc b/lib/maths/CPoissonMeanConjugate.cc
index bf897a6d85..3a45a9b62a 100644
--- a/lib/maths/CPoissonMeanConjugate.cc
+++ b/lib/maths/CPoissonMeanConjugate.cc
@@ -25,7 +25,6 @@
 #include <boost/math/distributions/negative_binomial.hpp>
 #include <boost/math/distributions/normal.hpp>
 #include <boost/math/distributions/poisson.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 
 #include <algorithm>
@@ -522,48 +521,43 @@ CPoissonMeanConjugate::jointLogMarginalLikelihood(const TDouble1Vec& samples,
     //   a is the prior gamma shape
     //   b is the prior gamma rate
 
-    try {
-        // Calculate the statistics we need for the calculation.
-        double numberSamples = 0.0;
-        double sampleSum = 0.0;
-        double sampleLogFactorialSum = 0.0;
-
-        for (std::size_t i = 0u; i < samples.size(); ++i) {
-            double n = maths_t::countForUpdate(weights[i]);
-            double x = samples[i] + m_Offset;
-            if (x < 0.0) {
-                // Technically, the marginal likelihood is zero here
-                // so the log would be infinite. We use minus max
-                // double because log(0) = HUGE_VALUE, which causes
-                // problems for Windows. Calling code is notified
-                // when the calculation overflows and should avoid
-                // taking the exponential since this will underflow
-                // and pollute the floating point environment. This
-                // may cause issues for some library function
-                // implementations (see fe*exceptflag for more details).
-                result = boost::numeric::bounds<double>::lowest();
-                return maths_t::E_FpOverflowed;
-            }
+    // Calculate the statistics we need for the calculation.
+    double numberSamples = 0.0;
+    double sampleSum = 0.0;
+    double sampleLogFactorialSum = 0.0;
 
-            numberSamples += n;
-            sampleSum += n * x;
-            // Recall n! = Gamma(n + 1).
-            sampleLogFactorialSum += n * boost::math::lgamma(x + 1.0);
+    for (std::size_t i = 0u; i < samples.size(); ++i) {
+        double n = maths_t::countForUpdate(weights[i]);
+        double x = samples[i] + m_Offset;
+        if (x < 0.0) {
+            // Technically, the marginal likelihood is zero here
+            // so the log would be infinite. We use minus max
+            // double because log(0) = HUGE_VALUE, which causes
+            // problems for Windows. Calling code is notified
+            // when the calculation overflows and should avoid
+            // taking the exponential since this will underflow
+            // and pollute the floating point environment. This
+            // may cause issues for some library function
+            // implementations (see fe*exceptflag for more details).
+            result = boost::numeric::bounds<double>::lowest();
+            return maths_t::E_FpOverflowed;
         }
 
-        // Get the implied shape parameter for the gamma distribution
-        // including the samples.
-        double impliedShape = m_Shape + sampleSum;
-        double impliedRate = m_Rate + numberSamples;
-
-        result = boost::math::lgamma(impliedShape) + m_Shape * std::log(m_Rate) -
-                 impliedShape * std::log(impliedRate) - sampleLogFactorialSum -
-                 boost::math::lgamma(m_Shape);
-    } catch (const std::exception& e) {
-        LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
-        return maths_t::E_FpFailed;
+        numberSamples += n;
+        sampleSum += n * x;
+        // Recall n! = Gamma(n + 1).
+        sampleLogFactorialSum += n * std::lgamma(x + 1.0);
     }
 
+    // Get the implied shape parameter for the gamma distribution
+    // including the samples.
+    double impliedShape = m_Shape + sampleSum;
+    double impliedRate = m_Rate + numberSamples;
+
+    result = std::lgamma(impliedShape) + m_Shape * std::log(m_Rate) -
+             impliedShape * std::log(impliedRate) - sampleLogFactorialSum -
+             std::lgamma(m_Shape);
+
     maths_t::EFloatingPointErrorStatus status = CMathsFuncs::fpStatus(result);
     if (status & maths_t::E_FpFailed) {
         LOG_ERROR(<< "Failed to compute log likelihood");
diff --git a/lib/maths/CRadialBasisFunction.cc b/lib/maths/CRadialBasisFunction.cc
index f9fee6d8de..91b30e8d0c 100644
--- a/lib/maths/CRadialBasisFunction.cc
+++ b/lib/maths/CRadialBasisFunction.cc
@@ -11,7 +11,6 @@
 #include <maths/CTools.h>
 
 #include <boost/math/constants/constants.hpp>
-#include <boost/math/special_functions/erf.hpp>
 
 #include <cmath>
 
@@ -33,7 +32,7 @@ double gaussianSquareDerivative(double x, double centre, double scale) {
     double r = scale * (x - centre);
     return scale *
            (boost::math::double_constants::root_two_pi *
-                boost::math::erf(boost::math::double_constants::root_two * r) -
+                std::erf(boost::math::double_constants::root_two * r) -
             4.0 * r * std::exp(-2.0 * r * r)) /
            4.0;
 }
@@ -51,7 +50,7 @@ double gaussianProduct(double x, double centre1, double centre2, double scale1,
     double d = scale1 * scale2 * (centre2 - centre1);
 
     return boost::math::double_constants::root_pi * std::exp(-d * d / (scale * scale)) *
-           boost::math::erf(scale * (x - m)) / (2.0 * scale);
+           std::erf(scale * (x - m)) / (2.0 * scale);
 }
 
 //! The indefinite integral
@@ -144,8 +143,7 @@ double CGaussianBasisFunction::mean(double a, double b, double centre, double sc
     }
 
     return std::max(boost::math::double_constants::root_pi / 2.0 / scale *
-                        (boost::math::erf(scale * (b - centre)) -
-                         boost::math::erf(scale * (a - centre))) /
+                        (std::erf(scale * (b - centre)) - std::erf(scale * (a - centre))) /
                         (b - a),
                     0.0);
 }
diff --git a/lib/maths/CStatisticalTests.cc b/lib/maths/CStatisticalTests.cc
index 9540622ebc..c53511700a 100644
--- a/lib/maths/CStatisticalTests.cc
+++ b/lib/maths/CStatisticalTests.cc
@@ -22,6 +22,7 @@
 #include <boost/range.hpp>
 
 #include <algorithm>
+#include <cmath>
 
 namespace ml {
 namespace maths {
@@ -70,7 +71,7 @@ double CStatisticalTests::leftTailFTest(double x, double d1, double d2) {
     if (x < 0.0) {
         return 0.0;
     }
-    if (boost::math::isinf(x)) {
+    if (std::isinf(x)) {
         return 1.0;
     }
     try {
@@ -87,7 +88,7 @@ double CStatisticalTests::rightTailFTest(double x, double d1, double d2) {
     if (x < 0.0) {
         return 1.0;
     }
-    if (boost::math::isinf(x)) {
+    if (std::isinf(x)) {
         return 0.0;
     }
     try {
diff --git a/lib/maths/CTools.cc b/lib/maths/CTools.cc
index 364c941dc7..925c2f8c18 100644
--- a/lib/maths/CTools.cc
+++ b/lib/maths/CTools.cc
@@ -32,7 +32,6 @@
 #include <boost/math/distributions/poisson.hpp>
 #include <boost/math/distributions/students_t.hpp>
 #include <boost/math/special_functions/digamma.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/tools/precision.hpp>
 #include <boost/optional.hpp>
 
@@ -1373,8 +1372,8 @@ double CTools::SIntervalExpectation::operator()(const normal& normal_, double a,
     double b_ = b == POS_INF ? b : (b - mean) / s;
     double expa = a_ == NEG_INF ? 0.0 : std::exp(-a_ * a_);
     double expb = b_ == POS_INF ? 0.0 : std::exp(-b_ * b_);
-    double erfa = a_ == NEG_INF ? -1.0 : boost::math::erf(a_);
-    double erfb = b_ == POS_INF ? 1.0 : boost::math::erf(b_);
+    double erfa = a_ == NEG_INF ? -1.0 : std::erf(a_);
+    double erfb = b_ == POS_INF ? 1.0 : std::erf(b_);
 
     if (erfb - erfa < std::sqrt(EPSILON)) {
         return expa == expb ? (a + b) / 2.0 : (a * expa + b * expb) / (expa + expb);
@@ -1405,8 +1404,8 @@ operator()(const lognormal& logNormal, double a, double b) const {
     double s = std::sqrt(2.0) * scale;
     double a_ = loga == NEG_INF ? NEG_INF : (loga - location) / s;
     double b_ = logb == POS_INF ? POS_INF : (logb - location) / s;
-    double erfa = loga == NEG_INF ? -1.0 : boost::math::erf((loga - c) / s);
-    double erfb = logb == POS_INF ? 1.0 : boost::math::erf((logb - c) / s);
+    double erfa = loga == NEG_INF ? -1.0 : std::erf((loga - c) / s);
+    double erfb = logb == POS_INF ? 1.0 : std::erf((logb - c) / s);
 
     if (erfb - erfa < std::sqrt(EPSILON)) {
         double expa = loga == NEG_INF ? 0.0 : std::exp(-a_ * a_);
@@ -1415,8 +1414,8 @@ operator()(const lognormal& logNormal, double a, double b) const {
                             : (expa + expb) / (expa / a + expb / b);
     }
 
-    double erfa_ = a_ == NEG_INF ? -1.0 : boost::math::erf(a_);
-    double erfb_ = b_ == POS_INF ? 1.0 : boost::math::erf(b_);
+    double erfa_ = a_ == NEG_INF ? -1.0 : std::erf(a_);
+    double erfb_ = b_ == POS_INF ? 1.0 : std::erf(b_);
     return mean * (erfb - erfa) / (erfb_ - erfa_);
 }
 
@@ -1845,7 +1844,7 @@ double CTools::differentialEntropy(const gamma& gamma_) {
 
     double shape = gamma_.shape();
     double scale = gamma_.scale();
-    return shape + std::log(scale) + boost::math::lgamma(shape) +
+    return shape + std::log(scale) + std::lgamma(shape) +
            (1 - shape) * boost::math::digamma(shape);
 }
 
diff --git a/lib/maths/ProbabilityAggregators.cc b/lib/maths/ProbabilityAggregators.cc
index b4449f5ecb..d7c07afe25 100644
--- a/lib/maths/ProbabilityAggregators.cc
+++ b/lib/maths/ProbabilityAggregators.cc
@@ -17,10 +17,11 @@
 
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/distributions/normal.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 #include <boost/range.hpp>
 
+#include <cmath>
+
 namespace ml {
 namespace maths {
 namespace {
@@ -144,8 +145,8 @@ class CNumericalLogProbabilityOfMFromNExtremeSamples {
         double result;
         CLogIntegrand f(m_P, m_Corrections, m_N, m_P.size(), 1u);
         CIntegration::logGaussLegendre<CIntegration::OrderThree>(f, 0, m_P[0], result);
-        result += boost::math::lgamma(static_cast<double>(m_N) + 1.0) -
-                  boost::math::lgamma(static_cast<double>(m_N - m_P.size()) + 1.0);
+        result += std::lgamma(static_cast<double>(m_N) + 1.0) -
+                  std::lgamma(static_cast<double>(m_N - m_P.size()) + 1.0);
         return result;
     }
 
@@ -496,7 +497,7 @@ bool CLogJointProbabilityOfLessLikelySamples::calculateLowerBound(double& result
     try {
         double logx = std::log(x);
         double p = std::floor(s - 1.0);
-        double logPFactorial = boost::math::lgamma(p + 1.0);
+        double logPFactorial = std::lgamma(p + 1.0);
         double m = std::floor(std::min(x, p) + 0.5);
         double logm = std::log(m);
 
@@ -546,7 +547,7 @@ bool CLogJointProbabilityOfLessLikelySamples::calculateLowerBound(double& result
 
         double logSum = logPFactorial - p * logx + std::max(b1, b2);
 
-        bound = (s - 1.0) * logx - x + logSum - boost::math::lgamma(s);
+        bound = (s - 1.0) * logx - x + logSum - std::lgamma(s);
 
         LOG_TRACE(<< "s = " << s << ", x = " << x << ", p = " << p
                   << ", m = " << m << ", b1 = " << b1 << ", b2 = " << b2
@@ -657,11 +658,11 @@ bool CLogJointProbabilityOfLessLikelySamples::calculateUpperBound(double& result
             b1 = std::log(p + 1);
         }
 
-        double b2 = boost::math::lgamma(p + 1.0) - p * std::log(x) + x;
+        double b2 = std::lgamma(p + 1.0) - p * std::log(x) + x;
 
         double logSum = std::min(b1, b2);
 
-        bound = (s - 1.0) * std::log(x) - x + logSum - boost::math::lgamma(s);
+        bound = (s - 1.0) * std::log(x) - x + logSum - std::lgamma(s);
 
         LOG_TRACE(<< "s = " << s << ", x = " << x << ", b1 = " << b1 << ", b2 = " << b2
                   << ", log(sum) = " << logSum << ", bound = " << bound);
@@ -876,8 +877,8 @@ bool CLogProbabilityOfMFromNExtremeSamples::calculate(double& result) {
 
     if (M > 1) {
         double logScale = static_cast<double>(M) * std::log(2.0) +
-                          boost::math::lgamma(static_cast<double>(N + 1)) -
-                          boost::math::lgamma(static_cast<double>(N - M + 1)) + logLargestCoeff;
+                          std::lgamma(static_cast<double>(N + 1)) -
+                          std::lgamma(static_cast<double>(N - M + 1)) + logLargestCoeff;
         LOG_TRACE(<< "log(scale) = " << logScale);
 
         double sum = 0.0;
diff --git a/lib/maths/unittest/CIntegerToolsTest.cc b/lib/maths/unittest/CIntegerToolsTest.cc
index dc814f92a8..8a1a7a4ce7 100644
--- a/lib/maths/unittest/CIntegerToolsTest.cc
+++ b/lib/maths/unittest/CIntegerToolsTest.cc
@@ -13,7 +13,6 @@
 
 #include <test/CRandomNumbers.h>
 
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/numeric/conversion/bounds.hpp>
 #include <boost/range.hpp>
 
@@ -183,10 +182,9 @@ void CIntegerToolsTest::testBinomial() {
             LOG_DEBUG(<< "j = " << j << ", n = " << n[i]
                       << ", (n j) = " << maths::CIntegerTools::binomial(n[i], j));
 
-            double expected =
-                std::exp(boost::math::lgamma(static_cast<double>(n[i] + 1)) -
-                         boost::math::lgamma(static_cast<double>(n[i] - j + 1)) -
-                         boost::math::lgamma(static_cast<double>(j + 1)));
+            double expected = std::exp(std::lgamma(static_cast<double>(n[i] + 1)) -
+                                       std::lgamma(static_cast<double>(n[i] - j + 1)) -
+                                       std::lgamma(static_cast<double>(j + 1)));
             CPPUNIT_ASSERT_DOUBLES_EQUAL(
                 expected, maths::CIntegerTools::binomial(n[i], j), 1e-10);
         }
diff --git a/lib/maths/unittest/CProbabilityAggregatorsTest.cc b/lib/maths/unittest/CProbabilityAggregatorsTest.cc
index 6802122a35..d01e417681 100644
--- a/lib/maths/unittest/CProbabilityAggregatorsTest.cc
+++ b/lib/maths/unittest/CProbabilityAggregatorsTest.cc
@@ -19,6 +19,8 @@
 #include <boost/optional.hpp>
 #include <boost/range.hpp>
 
+#include <cmath>
+
 using namespace ml;
 using namespace maths;
 using namespace test;
@@ -125,8 +127,8 @@ class CExpectedLogProbabilityOfMFromNExtremeSamples {
         TDoubleVec p(m_P.begin(), m_P.end());
         CLogIntegrand f(p, m_N, p.size(), 1u);
         CIntegration::logGaussLegendre<CIntegration::OrderTen>(f, 0, p[0], result);
-        result += boost::math::lgamma(static_cast<double>(m_N) + 1.0) -
-                  boost::math::lgamma(static_cast<double>(m_N - p.size()) + 1.0);
+        result += std::lgamma(static_cast<double>(m_N) + 1.0) -
+                  std::lgamma(static_cast<double>(m_N - p.size()) + 1.0);
         return result;
     }
 
@@ -282,7 +284,7 @@ void CProbabilityAggregatorsTest::testLogJointProbabilityOfLessLikelySamples() {
         double s = jointProbability.numberSamples() / 2.0;
         double x = jointProbability.distance() / 2.0;
 
-        double logP = logUpperIncompleteGamma(s, x) - boost::math::lgamma(s);
+        double logP = logUpperIncompleteGamma(s, x) - std::lgamma(s);
         LOG_DEBUG(<< "log(p) = " << logP);
 
         double lowerBound, upperBound;
@@ -328,7 +330,7 @@ void CProbabilityAggregatorsTest::testLogJointProbabilityOfLessLikelySamples() {
                     double x = jointProbability.distance() / 2.0;
                     LOG_DEBUG(<< "s = " << s << ", x = " << x);
 
-                    double logP = logUpperIncompleteGamma(s, x) - boost::math::lgamma(s);
+                    double logP = logUpperIncompleteGamma(s, x) - std::lgamma(s);
                     LOG_DEBUG(<< "log(p) = " << logP);
 
                     double lowerBound, upperBound;
@@ -348,7 +350,7 @@ void CProbabilityAggregatorsTest::testLogJointProbabilityOfLessLikelySamples() {
                     double s = jointProbability.numberSamples() / 2.0;
                     double x = jointProbability.distance() / 2.0;
 
-                    double logP = logUpperIncompleteGamma(s, x) - boost::math::lgamma(s);
+                    double logP = logUpperIncompleteGamma(s, x) - std::lgamma(s);
 
                     double lowerBound, upperBound;
                     CPPUNIT_ASSERT(logJointProbability.calculateLowerBound(lowerBound));
diff --git a/lib/maths/unittest/CSamplingTest.cc b/lib/maths/unittest/CSamplingTest.cc
index 4a30353cd0..c00bbfa231 100644
--- a/lib/maths/unittest/CSamplingTest.cc
+++ b/lib/maths/unittest/CSamplingTest.cc
@@ -12,9 +12,9 @@
 #include <maths/CBasicStatistics.h>
 #include <maths/CSampling.h>
 
-#include <boost/math/special_functions/gamma.hpp>
 #include <boost/range.hpp>
 
+#include <cmath>
 #include <numeric>
 
 using TDoubleVec = std::vector<double>;
@@ -28,11 +28,11 @@ using TDoubleVecVec = std::vector<TDoubleVec>;
 
 double multinomialProbability(const TDoubleVec& probabilities, const TSizeVec& counts) {
     std::size_t n = std::accumulate(counts.begin(), counts.end(), std::size_t(0));
-    double logP = boost::math::lgamma(static_cast<double>(n + 1));
+    double logP = std::lgamma(static_cast<double>(n + 1));
     for (std::size_t i = 0u; i < counts.size(); ++i) {
         double ni = static_cast<double>(counts[i]);
         if (ni > 0.0) {
-            logP += ni * std::log(probabilities[i]) - boost::math::lgamma(ni + 1.0);
+            logP += ni * std::log(probabilities[i]) - std::lgamma(ni + 1.0);
         }
     }
     return std::exp(logP);
diff --git a/lib/maths/unittest/CToolsTest.cc b/lib/maths/unittest/CToolsTest.cc
index 1cec2f3000..e7cef00356 100644
--- a/lib/maths/unittest/CToolsTest.cc
+++ b/lib/maths/unittest/CToolsTest.cc
@@ -237,8 +237,8 @@ double numericalProbabilityOfLessLikelySample(const boost::math::beta_distributi
 
     double xmin = 1000.0 * std::numeric_limits<double>::min();
     if (a >= 1.0 && fx < CTools::safePdf(beta, xmin)) {
-        return std::exp(a * std::log(xmin) - std::log(a) + boost::math::lgamma(a + b) -
-                        boost::math::lgamma(a) - boost::math::lgamma(b)) +
+        return std::exp(a * std::log(xmin) - std::log(a) + std::lgamma(a + b) -
+                        std::lgamma(a) - std::lgamma(b)) +
                CTools::safeCdfComplement(beta, x);
     }
 
@@ -1198,19 +1198,20 @@ void CToolsTest::testLgamma() {
     }
 
     double result;
-    CPPUNIT_ASSERT(maths::CTools::lgamma(0, result));
+    CPPUNIT_ASSERT(maths::CTools::lgamma(0, result, false));
     CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
 
     CPPUNIT_ASSERT((maths::CTools::lgamma(0, result, true) == false));
     CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
 
-    CPPUNIT_ASSERT((maths::CTools::lgamma(-1, result)));
+    CPPUNIT_ASSERT((maths::CTools::lgamma(-1, result, false)));
     CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
 
     CPPUNIT_ASSERT((maths::CTools::lgamma(-1, result, true) == false));
     CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
 
-    CPPUNIT_ASSERT((maths::CTools::lgamma(std::numeric_limits<double>::max() - 1, result)));
+    CPPUNIT_ASSERT((maths::CTools::lgamma(std::numeric_limits<double>::max() - 1,
+                                          result, false)));
     CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
 
     CPPUNIT_ASSERT((maths::CTools::lgamma(std::numeric_limits<double>::max() - 1,