//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [x] = CL_quantile(X,p)
// Quantiles of a data set (inverse of cumulative distribution function) 
//
// Calling Sequence
// x = CL_quantile(X,p)
//
// Description
// <itemizedlist><listitem>
// <p>Evaluates quantiles of a data set.</p>
// <p>Let F be the empirical cumulative distribution function associated with the data set <b>X</b>. 
// The quantity (x) that is computed is such that: </p>
// <p>F(x) = p (p: probability in [0,1]). </p>
// <p></p>
// <p>The method used to calculate sample quantiles is the same as in R-project (method number 5):</p>
// <p>For a data set containing N elements, quantiles are computed as follows:</p>
// <p>1) The sorted values in X are considered as the (0.5/N), (1.5/N),
//       ..., ((N-0.5)/N) quantiles.</p>
// <p>2) Linear interpolation is used to compute quantiles for
//       probabilities between (0.5/N) and ((N-0.5)/N)</p>
// <p>3) The minimum or maximum values in X are assigned to quantiles
//       for probabilities outside that range.</p>
// <p></p>
// <p>Note: This method is used in Matlab's quantile function as well.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// X: Data set (1xN)
// p: Probabilities (1xP)
// x: Quantiles (1xP)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) R: A Language and Environment for Statistical Computing; http://cran.r-project.org/doc/manuals/fullrefman.pdf. 
//
// Examples
// // Random gaussian sample
// N = 10000;
// m = 0;  // mean
// sd = 1; // standard deviation
// X = grand(1,N,"nor",m,sd);
//
// // Probability values
// p = [0.6827, 0.9545, 0.9973];
// x = CL_quantile(abs(X),p); // close to [1,2,3]*sd
// 
    
if (find(p < 0 | p > 1) <> [])
  error("Invalid probability values (must be in [0,1])"); 
end

if (size(X,1) <> 1)
  error("Invalid size for X (must be a row vector)");
end

X = gsort(X,"c","i");
N = length(X);

// (ptab,X) such that F is a piecewise linear function where 
// the knots are the values midway through the steps of the empirical cdf
ptab = [0, (0.5+(0:N-1))/N, 1];
X = [X(1), X, X($)];
 
// Linear interpolation 
x = CL_interpLin(ptab, X, p);
  
endfunction