diff --git a/include/evalhyd/detail/masks.hpp b/include/evalhyd/detail/masks.hpp
index 7a818ba5d76f605b60ecedd42d73ed91969891e6..60b94f229a8d2ca94045e90999371c7b3abc3486 100644
--- a/include/evalhyd/detail/masks.hpp
+++ b/include/evalhyd/detail/masks.hpp
@@ -19,8 +19,6 @@
 #include <xtensor/xsort.hpp>
 #include <xtensor/xindex_view.hpp>
 
-#include "maths.hpp"
-
 
 typedef std::map<std::string, std::vector<std::vector<std::string>>> msk_tree;
 
@@ -234,6 +232,8 @@ namespace evalhyd
                     auto q = get_q();
 
                     // define lambda function to precompute mean/median/quantile
+
+
                     auto get_val =
                             [&](const std::string& str, const std::string& num)
                     {
@@ -251,9 +251,8 @@ namespace evalhyd
                         }
                         else  // (str == "quantile")
                         {
-                            return evalhyd::maths::quantile(q, std::stod(num));
+                            return xt::quantile(q, {std::stod(num)})();
                         }
-
                     };
 
                     // preprocess conditions to identify special cases
diff --git a/include/evalhyd/detail/maths.hpp b/include/evalhyd/detail/maths.hpp
index be7db6af26fbe0bb2ffcacd403b4b32b7de4700a..d1601b736454ac85d9d26188fd6742967b1b8581 100644
--- a/include/evalhyd/detail/maths.hpp
+++ b/include/evalhyd/detail/maths.hpp
@@ -27,54 +27,6 @@ namespace evalhyd
                     xt::nanmean(xt::square(xt::abs(arr - mean_arr)), -1)
             );
         }
-
-        // function to calculate quantile on 1-dim expressions
-        inline double quantile(const xt::xtensor<double, 1>& t, double p)
-        {
-            std::size_t n = t.size();
-
-            // compute virtual index
-            auto virtual_idx = (double(n) - 1) * p;
-
-            if (std::fmod(virtual_idx, 1) == 0)
-            // virtual index is an integer
-            {
-                std::size_t i = {static_cast<unsigned long long>(virtual_idx)};
-                std::array<std::size_t, 1> kth = {i};
-                auto values = xt::partition(t, kth);
-                return xt::mean(xt::view(values, xt::range(i, i + 1)))();
-            }
-            else
-            // virtual index is not an integer
-            {
-                // determine indices before and after virtual index
-                std::size_t prv_idx = std::floor(virtual_idx);
-                std::size_t nxt_idx = prv_idx + 1;
-
-                // deal with edge cases
-                if (virtual_idx > double(n) - 1)
-                {
-                    prv_idx = n - 1;
-                    nxt_idx = n - 1;
-                }
-                if (virtual_idx < 0)
-                {
-                    prv_idx = 0;
-                    nxt_idx = 0;
-                }
-
-                std::array<std::size_t, 2> kth = {prv_idx, nxt_idx};
-                auto values = xt::partition(t, kth);
-                auto vw = xt::view(values, xt::range(prv_idx, nxt_idx + 1));
-
-                // perform linear interpolation to determine quantile
-                auto gamma = virtual_idx - double(prv_idx);
-                auto prv = xt::amin(vw);
-                auto nxt = xt::amax(vw);
-
-                return (prv + (nxt - prv) * gamma)();
-            }
-        }
     }
 }
 
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index 9ccf2add4b15b937e84bade912c340ab07e4e0e5..d1f451eae73d7ff7b19b62eeaebd8d3b48b93e7a 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -10,8 +10,6 @@
 #include <xtensor/xindex_view.hpp>
 #include <xtensor/xsort.hpp>
 
-#include "../maths.hpp"
-
 
 namespace evalhyd
 {
@@ -61,31 +59,14 @@ namespace evalhyd
                 xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
 
                 // compute predictive interval bounds from quantiles
-
-                // TODO: replace with `xt::quantile` when available
-                for (std::size_t s = 0; s < n_sit; s++)
+                for (std::size_t i = 0; i < n_itv; i++)
                 {
-                    for (std::size_t l = 0; l < n_ldt; l++)
-                    {
-                        for (std::size_t i = 0; i < n_itv; i++)
-                        {
-                            for (std::size_t t = 0; t < n_tim; t++)
-                            {
-                                // lower bound
-                                xt::view(itv_bnds, s, l, i, 0, t) =
-                                        evalhyd::maths::quantile(
-                                                xt::view(q_prd, s, l, xt::all(), t),
-                                                quantiles(i, 0)
-                                        );
-                                // upper bound
-                                xt::view(itv_bnds, s, l, i, 1, t) =
-                                        evalhyd::maths::quantile(
-                                                xt::view(q_prd, s, l, xt::all(), t),
-                                                quantiles(i, 1)
-                                        );
-                            }
-                        }
-                    }
+                    auto q =  xt::quantile(q_prd, xt::view(quantiles, i), 2);
+
+                    xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
+                            xt::view(q, 0);
+                    xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
+                            xt::view(q, 1);
                 }
                 
                 return itv_bnds;
@@ -214,7 +195,7 @@ namespace evalhyd
 
                         // compute "climatology" interval
 
-                        // TODO: replace with `xt::quantile` when available
+                        // TODO: replace with `xt::nanquantile` when available
                         for (std::size_t s = 0; s < n_sit; s++)
                         {
                             for (std::size_t l = 0; l < n_ldt; l++)
@@ -228,14 +209,14 @@ namespace evalhyd
 
                                     // lower bound
                                     xt::view(clim_bnds, s, l, m, e, i, 0) =
-                                            evalhyd::maths::quantile(
-                                                    obs_filtered, quantiles(i, 0)
-                                            );
+                                            xt::quantile(
+                                                    obs_filtered, {quantiles(i, 0)}
+                                            )();
                                     // upper bound
                                     xt::view(clim_bnds, s, l, m, e, i, 1) =
-                                            evalhyd::maths::quantile(
-                                                    obs_filtered, quantiles(i, 1)
-                                            );
+                                            xt::quantile(
+                                                    obs_filtered, {quantiles(i, 1)}
+                                            )();
                                 }
                             }
                         }