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

#import <lib/densest>

//
// The next program gives examples of the use of EstDensity, in
// comparison with DrawDensity. The results are usually very similar,
// with the exception of Cauchy	data, because the default window
// width 1.06*sigma*n^-0.2 doesn't work there.
//


main()
{
	decl x, ret, t, fhat, f, fOx;
	decl time = timer();

    print("Processing N(0,1) ... ");
    x = rann(50000,1);
	ret = DensEst(x,0,0,0,0);
	t = ret[0][];
	fhat = ret[2][];
	f = densn(t);
	DrawXMatrix(0,fhat,{"DensEst"},t,0);
	DrawDensity(1,x',"DrawDensity",1,0,0);
	DrawXMatrix(2,f,{"N(0,1) density"},t,0);
	SetDrawWindow("N(0,1) data");

	ret = DensEst(x,0,4,0,0);
	t = ret[0][];
	fhat = ret[2][];
	f = densn(t);
	DrawXMatrix(3,fhat,{"DensEst"},t,0);
	DrawDensity(4,x',"DrawDensity",1,0,0);
	DrawXMatrix(5,f,{"N(0,1) density"},t,0);
	SetDrawWindow("N(0,1) data");
	ShowDrawWindow();

    print("U[0,1] ... ");
    x = ranu(50000,1);
	ret = DensEst(x,0,0,0,0);
	t = ret[0][];
	fhat = ret[2][];
	f = (t.>0).*(t.<1);
	DrawXMatrix(0,fhat,{"DensEst"},t,0);
	DrawDensity(1,x',"DrawDensity",1,0,0);
	DrawXMatrix(2,f,{"U[0,1] density"},t,0);
	DrawXMatrix(3,f,{"U[0,1] density"},t,0);
	SetDrawWindow("U[0,1] data");
	ShowDrawWindow();

    print("Chi^2(5) ... ");
	x = ranchi(50000,1,5);
	ret = DensEst(x,0,0,0,0);
	t = ret[0][];
	fhat = ret[2][];
	f = denschi(t,5);
	DrawXMatrix(0,fhat,{"DensEst"},t,0);
	DrawDensity(1,x',"DrawDensity",1,0,0);
	DrawXMatrix(2,f,{"Chi^2(5) density"},t,0);
	DrawXMatrix(3,f,{"Chi^2(5) density"},t,0);
	SetDrawWindow("Chi^2(5) data");
	ShowDrawWindow();

    print("Cauchy ... ");
    x = rant(50000,1,1);
	ret = DensEst(x,-10,10,0.12,0);      // note that default window width = infinity 
	t = ret[0][];
	fhat = ret[2][];
	f = 1.0./((t.^2+1)*M_PI);
	fOx = denst(t,1);
    DrawXMatrix(0,fhat,{"DensEst"},t,0);
    DrawDensity(1,x',"DrawDensity",1,0,0);
	DrawXMatrix(2,f,{"Cauchy density"},t,0);
	DrawXMatrix(3,fOx,{"Ox t(1) density"},t,0);
	SetDrawWindow("Cauchy data");
	ShowDrawWindow();

	println("time = ", timespan(time));
}

