diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index 9343093faebfb4cab62bb6fae3f2e4f9d075a02c..65c133d33d6577dcca92834ab2aa05d8030a103c 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -20,18 +20,16 @@ TEST(ProbabilistTests, TestBrier) {
     xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
     ifs.close();
 
-    // add "artificial" site dimension (size 1) to observations
-    auto obs = xt::view(observed, xt::newaxis(), xt::all());
-    // add "artificial" site dimension (size 1) and lead time dimension (size 1)
-    // to predictions
-    auto prd = xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all());
-
     // compute scores
     xt::xtensor<double, 1> thresholds = {690, 534, 445};
 
     std::vector<xt::xarray<double>> metrics =
             evalhyd::evalp(
-                    obs, prd, {"BS", "BSS", "BS_CRD", "BS_LBD"},
+                    // shape: (sites [1], time [t])
+                    xt::view(observed, xt::newaxis(), xt::all()),
+                    // shape: (sites [1], lead times [1], members [m], time [t])
+                    xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
+                    {"BS", "BSS", "BS_CRD", "BS_LBD"},
                     thresholds
             );
 
@@ -61,6 +59,53 @@ TEST(ProbabilistTests, TestBrier) {
     EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd));
 }
 
+TEST(ProbabilistTests, TestQuantiles) {
+    // read in data
+    std::ifstream ifs;
+    ifs.open("./data/q_obs.csv");
+    xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<int>(ifs));
+    ifs.close();
+
+    ifs.open("./data/q_prd.csv");
+    xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
+    ifs.close();
+
+    // compute scores
+    xt::xtensor<double, 1> thresholds = {690, 534, 445};
+
+    std::vector<xt::xarray<double>> metrics =
+            evalhyd::evalp(
+                    // shape: (sites [1], time [t])
+                    xt::view(observed, xt::newaxis(), xt::all()),
+                    // shape: (sites [1], lead times [1], members [m], time [t])
+                    xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
+                    {"QS", "CRPS"},
+                    thresholds
+            );
+
+    // check results
+    // Quantile scores
+    xt::xtensor<double, 4> qs =
+            {{{{ 345.91578 ,  345.069256,  343.129359,  340.709869,  338.281598,
+                 335.973535,  333.555157,  330.332426,  327.333539,  324.325996,
+                 321.190082,  318.175117,  315.122186,  311.97205 ,  308.644942,
+                 305.612169,  302.169552,  298.445956,  294.974648,  291.273807,
+                 287.724586,  284.101905,  280.235592,  276.21865 ,  272.501484,
+                 268.652733,  264.740168,  260.8558  ,  256.90329 ,  252.926292,
+                 248.931239,  244.986396,  240.662998,  236.328964,  232.089785,
+                 227.387089,  222.976008,  218.699975,  214.099678,  209.67252 ,
+                 205.189587,  200.395746,  195.2372  ,  190.080139,  185.384244,
+                 180.617858,  174.58323 ,  169.154093,  163.110932,  156.274796,
+                 147.575315}}}};
+    EXPECT_TRUE(xt::allclose(metrics[0], qs));
+
+    // Continuous ranked probability scores
+    xt::xtensor<double, 3> crps =
+            {{{ 257.412129}}};
+    EXPECT_TRUE(xt::allclose(metrics[1], crps));
+
+}
+
 TEST(ProbabilistTests, TestMasks) {
     // read in data
     std::ifstream ifs;
@@ -79,21 +124,27 @@ TEST(ProbabilistTests, TestMasks) {
 
     // compute scores using masks to subset whole record
     xt::xtensor<double, 1> thresholds = {690, 534, 445};
-    std::vector<std::string> metrics = {"BS", "BSS", "BS_CRD", "BS_LBD"};
+    std::vector<std::string> metrics =
+            {"BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS"};
 
     std::vector<xt::xarray<double>> metrics_masked =
             evalhyd::evalp(
+                    // shape: (sites [1], time [t])
                     xt::view(observed, xt::newaxis(), xt::all()),
+                    // shape: (sites [1], lead times [1], members [m], time [t])
                     xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
                     metrics,
                     thresholds,
+                    // shape: (subsets [1], time [t])
                     masks
             );
 
     // compute scores on pre-computed subset of whole record
     std::vector<xt::xarray<double>> metrics_subset =
             evalhyd::evalp(
+                    // shape: (sites [1], time [t-20])
                     xt::view(observed, xt::newaxis(), xt::range(20, _)),
+                    // shape: (sites [1], lead times [1], members [m], time [t-20])
                     xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _)),
                     metrics,
                     thresholds