/*--------------------------------------------------------------------------
 * 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 g_M = 128 points.
 *                                                                            
 * Syntax:
 *
 * DensEst(x, minx, maxx, h)                                          
 *     x     in: nx1 vector of observations                               
 *     minx  in: minimum x-value at which f is estimated;
 *               default (minx>=maxx): minx = min(x) - 3*h;             
 *     maxx  in: maximum x-value at which f is estimated
 *      		 default (minx>=maxx): maxx = max(x) + 3*h;
 *     h     in: window width
 *               default (h<=0): h = 1.06*sigma*n^-0.2
 *
 * 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 = g_M = 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 g_M = 128;

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

	cx = rows(vX);
	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;
}
