Index: tol/bmath/bstat/stat.cpp
===================================================================
--- tol/bmath/bstat/stat.cpp	(revision 678076951fff3e4eb0155adbc1338577af44dac5)
+++ tol/bmath/bstat/stat.cpp	(revision a106e90dc95a85e4d0b06fae1284069e93887d41)
@@ -383,4 +383,52 @@
     }
     return(result);
+}
+
+//--------------------------------------------------------------------
+void Quantile(const BArray<BDat>& vec, 
+              const BArray<BDat>& r, 
+                    BArray<BDat>& Qu)
+
+/*! Returns a set of q-Quantile statistics of a vector of real numbers
+ */
+//--------------------------------------------------------------------
+{
+  Qu.AllocBuffer(0);
+  if(!vec.Size()) { return; }
+  BArray<BDat> sorted(vec.Size());
+  int i, j;
+  for(i=j=0; i<vec.Size(); i++)
+  {
+    if(vec(i).IsKnown())
+    {
+      sorted(j++) = vec(i);
+    }
+  }
+  sorted.ReallocBuffer(j);
+  sorted.Sort(DatCmp);
+  Qu.AllocBuffer(r.Size());
+  int k;
+  for(k=0; k<r.Size(); k++)
+  {
+    BDat q = r[k];
+    if(q>1)  { q = 1; }
+    if(q<0)  { q = 0; }
+    BDat p = 1-q;
+    BDat Q = (1.0/q).Value();
+    BDat n = BDat(sorted.Size()-1)/BDat(Q);
+    BDat result;
+    if(n.IsInteger())
+    {
+      Qu[k] = sorted(int(n.Value()));
+    }
+    else
+    {
+      BInt nq = (BInt)Ceil (n).Value();
+      BInt np = (BInt)Floor(n).Value();
+      BDat sq = sorted(nq);
+      BDat sp = sorted(np);
+      Qu[k] = (q*sq)+(p*sp);
+    }
+  }
 }
 
Index: tol/bmath/bstat/tol_bstat.h
===================================================================
--- tol/bmath/bstat/tol_bstat.h	(revision 678076951fff3e4eb0155adbc1338577af44dac5)
+++ tol/bmath/bstat/tol_bstat.h	(revision a106e90dc95a85e4d0b06fae1284069e93887d41)
@@ -33,57 +33,61 @@
 		      const BArray<BDat>& vec);
 
-TOL_API BDat Sum	         (const BArray<BDat>& vec);
-TOL_API BDat Product	         (const BArray<BDat>& vec);
-TOL_API BDat Average	         (const BArray<BDat>& vec);
-TOL_API BDat HarmonicAverage     (const BArray<BDat>& vec);
-TOL_API BDat GeometricAverage    (const BArray<BDat>& vec);
-TOL_API BDat Moment	         (const BArray<BDat>& vec, BInt m);
-TOL_API BDat CenterMoment        (const BArray<BDat>& vec, BInt m);
-TOL_API BDat Varianze	         (const BArray<BDat>& vec);
-TOL_API BDat StandardDeviation   (const BArray<BDat>& vec);
+TOL_API BDat Sum	                (const BArray<BDat>& vec);
+TOL_API BDat Product	            (const BArray<BDat>& vec);
+TOL_API BDat Average	            (const BArray<BDat>& vec);
+TOL_API BDat HarmonicAverage      (const BArray<BDat>& vec);
+TOL_API BDat GeometricAverage     (const BArray<BDat>& vec);
+TOL_API BDat Moment	              (const BArray<BDat>& vec, BInt m);
+TOL_API BDat CenterMoment         (const BArray<BDat>& vec, BInt m);
+TOL_API BDat Varianze	            (const BArray<BDat>& vec);
+TOL_API BDat StandardDeviation    (const BArray<BDat>& vec);
 TOL_API BDat AsymmetryCoefficient (const BArray<BDat>& vec);
-TOL_API BDat Kurtosis	         (const BArray<BDat>& vec);
-TOL_API BDat Maximum	         (const BArray<BDat>& vec);
-TOL_API BDat Minimum	         (const BArray<BDat>& vec);
-TOL_API BDat Quantile	         (const BArray<BDat>& vec, BDat q);
-TOL_API BDat Median	         (const BArray<BDat>& vec);
-
-TOL_API BDat Covarianze (const BArray<BDat>& v1,
-		                 const BArray<BDat>& v2);
-TOL_API BDat Correlation (const BArray<BDat>& v1,
-		                  const BArray<BDat>& v2);
-
-TOL_API void Covarianze (const BMatrix<BDat>&    data,
-		                 BSymMatrix<BDat>& cor);
-TOL_API void Correlation (const BMatrix<BDat>&    data,
-		                  BSymMatrix<BDat>& cor);
-TOL_API void PartialCorrelation (const BMatrix<BDat>&    data,
-		                     	 BSymMatrix<BDat>& pcor);
-TOL_API void CorrelationF (const BSymMatrix<BDat>& cov,
-	                       BSymMatrix<BDat>& cor);
-TOL_API void PartialCorrelationF (const BSymMatrix<BDat>& cor,
-			                      BSymMatrix<BDat>& pcor);
-TOL_API void RegressionCoef (const BSymMatrix<BDat>& cov,
-		                     const BSymMatrix<BDat>& pcor,
-		                     BMatrix<BDat>&    regr);
-
-TOL_API BDat PearsonCoefficient (const BArray<BDat>& v1,
-		 	                    const BArray<BDat>& v2);
-
-TOL_API BDat R2Coefficient (const BArray<BDat>& v1, const BArray<BDat>& v2);
-
-TOL_API void AutoCov (const BArray<BDat>& vec, BArray<BDat>& cov, BInt order, 
-		      BInt demean=0);
-
-TOL_API void AutoCor (const BArray<BDat>& vec, BArray<BDat>& cor, BInt order, 
-		      BInt demean=0);
-
-TOL_API void DurbinAutoReg (const BArray<BDat>&  r, BMatrix<BDat>& FI);
-
-TOL_API void PartAutoCor   (const BArray<BDat>& cor, BArray<BDat>& p);
-TOL_API void InverseAutoCor(const BArray<BDat>& cor, BArray<BDat>& p);
-TOL_API void InvPartAutoCor(const BArray<BDat>& cor, BArray<BDat>& p);
-TOL_API BDat BoxPierceLjung(const BArray<BDat>& vec,
-		       	            BInt m, BInt p, BInt q);
+TOL_API BDat Kurtosis	            (const BArray<BDat>& vec);
+TOL_API BDat Maximum	            (const BArray<BDat>& vec);
+TOL_API BDat Minimum	            (const BArray<BDat>& vec);
+TOL_API BDat Quantile	            (const BArray<BDat>& vec, BDat q);
+TOL_API void Quantile             (const BArray<BDat>& vec, 
+                                   const BArray<BDat>& r, 
+                                         BArray<BDat>& Qu);
+TOL_API BDat Median	              (const BArray<BDat>& vec);
+TOL_API BDat Covarianze           (const BArray<BDat>& v1,
+		                               const BArray<BDat>& v2);
+TOL_API BDat Correlation          (const BArray<BDat>& v1,
+		                               const BArray<BDat>& v2);
+TOL_API void Covarianze           (const BMatrix<BDat>&    data,
+		                                     BSymMatrix<BDat>& cor);
+TOL_API void Correlation          (const BMatrix<BDat>&    data,
+		                                     BSymMatrix<BDat>& cor);
+TOL_API void PartialCorrelation   (const BMatrix<BDat>&    data,
+		                     	               BSymMatrix<BDat>& pcor);
+TOL_API void CorrelationF         (const BSymMatrix<BDat>& cov,
+	                                       BSymMatrix<BDat>& cor);
+TOL_API void PartialCorrelationF  (const BSymMatrix<BDat>& cor,
+			                                   BSymMatrix<BDat>& pcor);
+TOL_API void RegressionCoef       (const BSymMatrix<BDat>& cov,
+		                               const BSymMatrix<BDat>& pcor,
+		                                     BMatrix<BDat>&    regr);
+TOL_API BDat PearsonCoefficient   (const BArray<BDat>& v1,
+		 	                             const BArray<BDat>& v2);
+TOL_API BDat R2Coefficient        (const BArray<BDat>& v1, 
+                                   const BArray<BDat>& v2);
+TOL_API void AutoCov              (const BArray<BDat>& vec, 
+                                         BArray<BDat>& cov, 
+                                         BInt order, 
+		                                     BInt demean=0);
+TOL_API void AutoCor              (const BArray<BDat>& vec, 
+                                         BArray<BDat>& cor, 
+                                         BInt order, 
+		                                     BInt demean=0);
+TOL_API void DurbinAutoReg        (const BArray<BDat>&  r, 
+                                         BMatrix<BDat>& FI);
+TOL_API void PartAutoCor          (const BArray<BDat>& cor, 
+                                         BArray<BDat>& p);
+TOL_API void InverseAutoCor       (const BArray<BDat>& cor, 
+                                         BArray<BDat>& p);
+TOL_API void InvPartAutoCor       (const BArray<BDat>& cor, 
+                                         BArray<BDat>& p);
+TOL_API BDat BoxPierceLjung       (const BArray<BDat>& vec,
+		       	                             BInt m, BInt p, BInt q);
 
 //! To obtine the exponent of the best Box-Cox exponential transformation of a