diff --git a/QuantileRegression.m b/QuantileRegression.m
index b015d6f1..4431b0ed 100644
--- a/QuantileRegression.m
+++ b/QuantileRegression.m
@@ -308,15 +308,19 @@ I experimented with using DualLinearProgramming (provided by Mathematica) and wi
 
 Clear[NURBSBasis];
 
-Options[NURBSBasis] = {SplineClosed -> False};
+Options[NURBSBasis] = {SplineClosed -> False, "RelativeMargin" -> 0.05};
 
 NURBSBasis[data_?MatrixQ, n_Integer, opts : OptionsPattern[]] :=
     NURBSBasis[data, Table[n, Length @ data[[1]]] ];
 
 NURBSBasis[data_?MatrixQ, nsArg : { _Integer .. }, opts : OptionsPattern[]] :=
-     Module[{ns = nsArg, dim = Dimensions[data][[2]],
-      lsMinMaxes, cpts0, inds, cpts, lsBasis},
+    Module[{ns = nsArg, dim = Dimensions[data][[2]],
+      lsMinMaxes, relMargin, cpts0, inds, cpts, lsBasis},
 
+      (* Process options *)
+      relMargin = OptionValue[NURBSBasis, "RelativeMargin"];
+
+      (* Extend number of points per side spec *)
       Which[
         dim < Length[ns],
         ns = Take[ns, dim],
@@ -325,8 +329,17 @@ I experimented with using DualLinearProgramming (provided by Mathematica) and wi
         ns = Take[Flatten[Table[ns, dim]], dim]
       ];
 
+      (* Min-max per column *)
       lsMinMaxes = MinMax /@ Transpose[data];
 
+      (* The application of margins is needed to prevent instabilities of the basis value computations. *)
+      lsMinMaxes =
+          Map[ {
+            #[[1]] - (#[[2]] - #[[1]]) * relMargin,
+            #[[2]] + (#[[2]] - #[[1]]) * relMargin
+          }&, lsMinMaxes];
+
+      (* Make the basis *)
       cpts0 = Outer[{0} &, Sequence @@ Map[Table[i, {i, 0, 1, 1 / (# - 1)}] &, ns]];
       inds = Outer[List, Sequence @@ Map[Range, ns]];
 
@@ -340,6 +353,7 @@ I experimented with using DualLinearProgramming (provided by Mathematica) and wi
         inds
       ];
 
+      (* Result *)
       lsBasis
     ];