/*--------------------------------------------------------------------------
 * DensEst.ox - code for density estimation
 *
 *       code (C) H. Peter Boswijk 1996
 *		 (JAD): minor changes; made to use discretize().
 *
 * The function DensEst calculates a Kernel estimate of the probability
 * density function f(x) of a random sample, collected in a vector x.
 * The function is based on the normal kernel, and the density is
 * evaluated at cM = 128 points.
 *                                                                            
 * Syntax:
 *
 * DensEst(const vX, dMinX, dMaxX, dH, cM)                                          
 *     vX     in: nx1 vector of observations                               
 *     dMinX  in: minimum x-value at which f is estimated;
 *                default (dMinX>=dMaxX): minx = min(x) - 3*h;             
 *     dMaxX  in: maximum x-value at which f is estimated
 *      		  default (dMinX>=dMaxX): maxx = max(x) + 3*h;
 *     dH     in: window width h
 *                default (h<=0): h = 1.06*sigma*n^-0.2
 *     cM     in: number of points (<=0: 128)
 *
 * Return value
 *     Returns the 3xs matrix t|hist|f, where t is a row vector of points
 *     where f is evaluated, and hist and f are the corresponding histogram
 *     and density values. When minx and maxx have their default value,
 *     s = 128; but when minx and maxx are prespecified, then s = 122,
 *     because the density is also evaluated at 3 points below minc and 3
 *     points above maxc, which are not returned.
 *--------------------------------------------------------------------------*/

#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>

const decl DENSEST_M = 128;

DensEst(const vX, dMinX, dMaxX, dH, cM)
{
    decl cx, deflt, step, t, hist, s, ffthist, arg, fac, f, ret, cadd = 6, cm1;

	cx = rows(vX);
	if (cM <= 0)
		cM = DENSEST_M;
	cm1 = cM - 1;

	if (dH <= 0)
        dH = 1.06 * sqrt(varc(vX)) * (cx ^ -0.2);
	dH = double(dH);

	if (dMinX >= dMaxX)
    {
	    deflt = 1;
		f = limits(vX);				// get min,max from data
	    dMinX = f[0][0] - 3 * dH;	// data_min - 3h
	    dMaxX = f[1][0] + 3 * dH;	// data_max - 3h
	}
	else						  	// over selected interval
	{
		deflt = 0;
	    step = (dMaxX - dMinX) / (cM - 2 * cadd);
	    dMinX = dMinX - cadd * step;// add bit to left and right
	    dMaxX = dMaxX + cadd * step;
	}
    step = (dMaxX - dMinX) / cm1;
									// create horizontal axis
	t = dMinX + step * range(0, cm1);
    	            // discretize data, discards data outside range
	hist = discretize(vX', dMinX, dMaxX, cM, 0);
	hist /= (cx * step);
	
    ffthist = fft(hist);			// do the FFT, over [0,2pi)
	s = range(0, cm1) * M_2PI / (dMaxX - dMinX);
    fac = exp(-0.5 * sqr(dH * s));	// normal kernel

	f = fft(fac .* ffthist, 2);		// inverse real FFT 
	f = f .* (f .> 0);				// set negative part to zero
    ret = t | hist | f;

	if (deflt == 0)	// remove points outside original interval
        ret = ret[][cadd : cM-cadd-1];

	return ret;
}
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------
 * BivDensEst.ox - code for bivariate density estimation
 *
 *       code (C) H. Peter Boswijk 2003
 *		 (JAD): minor changes
 *
 * The function BivDensEst calculates a kernel estimate of the
 * bivariate probability density function f(x1,x2) of a random sample,
 * collected in an nx2 matrix mX.
 * The function is based on the Gaussian kernel, after standardizing
 * using the sample covariance matrix. The range of x1 and x2 may be set
 * by the user; the default choice is determined by the sample minimum
 * and maximum of x1 and x2. The window width may also be set by the user;
 * the default here is the optimal value for Gaussian data, which for
 * this dimension and kernel is 0.96n^(-1/6), see Silverman (1986,
 * section 4.3.2). The density is evaluated at (cM+1) x (cM+1) points,
 * where cM = 30 by default. 
 *
 * By setting the dH argument to -1 or -2, the Gaussian kernel is
 * replaced by the Epanechnikov and "squared Epanechnikov" kernel
 * respectively, see Silverman (1986, p.76, K_e and K_2),
 * with corresponding optimal window widths. This is slightly faster.
 *                                                                            
 * Syntax:
 *
 * BivDensEst(mX, dMinX, dMaxX, dH, cM)
 *     mX     in: nx2 matrix of observations                               
 *     dMinX  in: 2x1 vector of minimum x-values at which f is estimated
 *                default (dMinX==dMaxX==0): dMinX = minc(mX);             
 *     dMaxX  in: 2x1 vector of maximum x-values at which f is estimated
 *      	 	  default (dMinX==dMaxX==0): dMaxX = maxc(mX);
 *     dH     in: window width
 *                default (dH== 0): Gaussian kernel, h = 0.96*n^(-1/6).
 *                default (dH==-1): Epanechnikov, h = 1.77*n^(-1/6).
 *                default (dH==-2): squared Epanechnikov, h = 2.04*n^(-1/6).
 *
 * Returns the array {t, f}, with
 *     t 	  : 2x(cM+1) matrix of points of x1 and x2 where f is evaluated;
 *     f      : (cM+1)x(cM+1) matrix constaining the corresponding density values.
 *
 *--------------------------------------------------------------------------*/

const decl DENSEST_M2 = 30;

BivDensEst(mX, dMinX, dMaxX, dH, cM)
{
    decl n, step, t, f, chol, fact, i, j, dx, xtx, logdet, sqrtdet, ti, mxtx;
	decl kernel = 0;	// 0: Gaussian (default), 1: Epanechnikov, 2: squared Epanechnikov

	n = rows(mX);
	if (columns(mX) != 2)
		return {<>, <>};

	if (cM <= 0)
		cM = DENSEST_M2;
		
	if (dH == -1)
	{	dH = 1.77 * n ^ (-1/6);			// for Epanechnikov kernel
		kernel = 1;
	}
	else if (dH == -2)
	{	dH = 2.04 * n ^ (-1/6);			// for squared Epanechnikov kernel
		kernel = 2;
	}
	else if (dH <= 0)
		dH = 0.96 * n ^ (-1/6);			// optimal h for Gaussian kernel

	if (dMinX == 0 && dMaxX == 0)
    {
	    dMinX = minc(mX)';
	    dMaxX = maxc(mX)';
	}
    step = (dMaxX - dMinX) / cM;
	t = dMinX + step * range(0, cM);

	chol = choleski(invertsym(variance(mX), &logdet));
	sqrtdet = exp(-logdet / 2);			// sqrt(determinant(inv(S)))
	fact = (1/M_2PI) * sqrtdet / sqr(dH);
	if (kernel == 1)
		fact *= 4;						// for Epanechnikov instead of Gaussian kernel
	else if (kernel == 2)
		fact *= 6;						// for squared Epanechnikov instead of Gaussian kernel
	
	chol /= dH;							// move window scaling outside loop
	mX *= chol;							// move standardization out
	f = zeros(cM + 1, cM + 1);			// destination for kernel estimates

	for (i = 0; i <= cM; ++i)
	{
		ti = t[0][i];
		if (kernel == 1)				// Epanechnikov kernel
		{
			for (j = 0; j <= cM; ++j)
			{
				xtx = sumsqrr((ti ~ t[1][j]) * chol - mX);
				f[i][j] = meanc((1-xtx).*(xtx.<1));
			}
		}
		else if (kernel == 2)			// squared Epanechnikov kernel
		{
			for (j = 0; j <= cM; ++j)
			{
				xtx = sumsqrr((ti ~ t[1][j]) * chol - mX);
				f[i][j] = meanc(sqr((1-xtx).*(xtx.<1)));
			}
		}
		else							// Gaussian kernel
		{
			for (j = 0; j <= cM; ++j)
			{
				xtx = sumsqrr((ti ~ t[1][j]) * chol - mX);
				f[i][j] = meanc(exp(-xtx/2));
			}
		}
	}
	f *= fact;

	return {t, f};
}
/*--------------------------------------------------------------------------*/
