/*--------------------------------------------------------------------------
 * Quantile.ox - simple code to compute quantiles
 *
 *       code (C) Jurgen Doornik 1994-1995
 *
 * This illustrates how to compute quantiles of a probability
 * distribution, in case the quantiles are not avalaible in Ox.
 *
 * Syntax:
 *
 *	Quantile(const p, const n1, const n2, const fnDens, const fnProb)
 *		p[cP]    in:  vector with probabilities
 *		n1		 in:  first degrees of freedom
 *		n2		 in:  second degrees of freedom
 *      fnDens   in:  function computing density, of form;
 *					  fnDens(const x, const n1, const n2);
 *      fnProb   in:  function computing cumulative density, of form;
 *					  fnProb(const x, const n1, const n2);
 *	Returns the quantiles
 *--------------------------------------------------------------------------*/

#include <oxstd.h>

const decl QUANT_MAXIT = 40;
const decl QUANT_EPS   = 1e-10;


AQuantile(const p, const n1, const n2, const fnDens, const fnProb)
{
    decl i, t = 0.5, f, df, dt;

    for (i = 0; i < QUANT_MAXIT; i++)
    {
        f = fnProb(t, n1, n2);
        df = fnDens(t, n1, n2);
        dt = (p - f) / df;
        t += dt;
        if ( fabs(p - f) < QUANT_EPS && fabs(dt) < QUANT_EPS)
            break;
    }

return t;
}
Quantile(const p, const n1, const n2, const fnDens, const fnProb)
{
    decl i, t = p;

    for (i = 0; i < sizerc(p); ++i)
        t[i] = AQuantile(p[i], n1, n2, fnDens, fnProb);

return t;
}

