// goodness-of-fit.sci
// code scilab (7 janvier 2011)

// exec kolmogorov.sci; // a appeler depuis le fichier central (pb windows)

// la fonction suivante calcule la statistique de Kolmogorov--Smirnov
// lorsque la fonction de r\'epartition cible est donn\'ee par $F$.
// syntaxe : KSstat(<\'echantillon>, <fonction de r\'epartition>)
// exemple :
// --> deff("y = F(x)", "y = x");
// --> k = KSstat(u, F);
// --> p = 1-kolmogorovcdf(k, n);

function ks = KSstat(u, F)
	// local v n i;
	n = size(u,1); v = -gsort(-u); ks = 0;
	for i = 1:n; ks = max(ks, F(v(i))-(i-1)/n, i/n-F(v(i))); end
endfunction

// x est un vecteur colonne contenant les observations
// F est la fonction de r\'epartition `cible'' donn\'ee
// par exemple sous la forme
//
// deff("y = F(x)", "y = min(max(x, 0), 1)");
//
// Exemple :
// N = 100; u = zeros(N, 1);
// for i = 1:N; u(i) = rand(); end
// deff("y = F(x)", "y = min(max(x, 0), 1)");
// cgoftest(u, F);

// Si x est une matrice de dimensions $n\times m$ avec $m>1$,
// alors un second niveau de test est effectu\'e : celui de
// l'uniformit\'e des $p$-valeurs correspondant aux $m$
// \'echantillons x(:, i), $1\leq i\leq m$.

function cgoftest(x, F)
	// local xs x_n x_m statistic pvaleur Fxs;
	mprintf("cgoftest :\n");
	x_n = size(x,1); x_m = size(x,2);
	//
	// Test de Kolmogorov--Smirnov
	//
	mprintf("Test de Kolmogorov-Smirnov\n");
	mprintf("Taille de l''echantillon n = %d\n", x_n);
	statistic = KSstat(x(:, 1), F);
	mprintf("Statistique de test = %f\n", statistic);
	pvaleur = 1-klmcdf(statistic, x_n);
	mprintf("p-valeur approchee = %f\n", pvaleur);
	if x_n <= 100 then pvaleur = 1-kolmogorovcdf(statistic, x_n);
		mprintf("p-valeur exacte = %f\n", pvaleur); end
	//
	// Graphiques
	//
	mprintf("cgoftest, graphiques...\n");
	clf();
	xs = -gsort(-x(:,1)); rang = [1:x_n]'/x_n;
	for i = 1:x_n; Fxs(i) = F(xs(i)); end
	plot2d(xs, Fxs); plot2d2(xs, rang);
	xtitle("Simulations et loi theorique");
	legend(["distribution theorique"; "distribution observee"], 4);
	//
	// Deuxi\`eme niveau de test
	//
	// local Statistic Pvaleur uniformcdf;
	if x_m > 1 then
	    halt("taper return pour continuer");
	    mprintf("Test d''uniformite des p-valeurs :"\n);
		statistic = []; pvaleur = [];
		for i = 1:x_m;
			statistic = [statistic; KSstat(x(:, i), F)];
			if x_n <= 100 then
				pvaleur = [pvaleur; 1-kolmogorovcdf(statistic(i), x_n)];
			else
				pvaleur = [pvaleur; 1-klmcdf(statistic(i), x_n)];
			end
		end
		deff("p = uniformcdf(x)", "y = min(1, max(0, x))");
		Statistic = KSstat(pvaleur, uniformcdf);
		mprintf("Statistique de test = %f\n", Statistic);
		Pvaleur = 1-klmcdf(Statistic, x_m);
		mprintf("p-valeur approchee = %f\n", Pvaleur);
		if x_m <= 100 then
			Pvaleur = 1-kolmogorovcdf(Statistic, x_m);
			mprintf("p-valeur exacte = %f\n", Pvaleur);
		end
		clf(); plot(pvaleur);
		xtitle("Simulations et loi theorique");
		legend(["suite des p-valeurs"]);
		// halt("taper return pour continuer");
	end 
	mprintf("cgoftest, fin.\n");
endfunction

// attention aux ";" sous windows

// x est un vecteur colonne contenant les observations
// xtheo est le vecteur colonne des valeurs possibles theoriques
// ptheo celui des probabilites correspondantes.

// Exemple :
// N = 1000; p = 0.25;
// for i = 1:N; if rand() > 1-p then x(i) = 1; else x(i) = 0; end end
// dgoftest(x,[0; 1],[1-p; p]);

function dgoftest(x, xtheo, ptheo)
	// local i j ptheo_n x_n pobs statistic pvaleur total;
	mprintf("dgoftest :\n");
	//
	// Comptabilit\'e
	//
	x_n = size(x,1); ptheo_n = size(ptheo,1); pobs = zeros(ptheo_n, 1);
	for i = 1:ptheo_n
		for j = 1:x_n
			if x(j) == xtheo(i) then pobs(i) = pobs(i)+1; end
		end
	end
	total = sum(pobs); pobs = pobs/x_n;
	//
	// Test du khi-deux
	//
	// noter le redoublement des apostrophes dans les cha\^\i nes suivantes
	mprintf("Test du khi-deux d''adequation\n");
	mprintf("Taille de l''echantillon = %d\n", x_n);
	mprintf("Nombre theorique de valeurs = %d\n", ptheo_n);
	if total < x_n then
		mprintf("Des valeurs observees sont hors des valeurs theoriques.\n");
		mprintf("L''hypothese d''adequation est trivialement rejetee.\n");
	else
		statistic = 0; flag = %f;// bool\'een, false
		for i = 1:ptheo_n;
			statistic = statistic+(pobs(i)-ptheo(i))^2/ptheo(i);
			if pobs(i)*x_n < 5 then flag = %t; end
		end
		statistic = x_n*statistic;
		mprintf("Statistique du khi-deux d''adequation = %f\n", statistic);
		if x_n < 30 | flag then // le ``ou'' logique est la barre verticale
			mprintf("Les conditions d''application du test asymptotique\n");
			mprintf("ne sont pas satisfaites. Cependant,\n");
		end
		pvaleur = 1-cdfchi("PQ", statistic, ptheo_n-1);
		mprintf("p-valeur asymtotique = %f\n", pvaleur);
	end
	//
	// Graphiques
	//
	printf("dgoftest, graphique et fin.\n");
	clf();
	bar(xtheo, [ptheo, pobs]);
	xtitle("Simulations et loi theorique");
	legend(["distribution theorique"; "distribution observee"]);
	a = gca(); a.grid = [-1,3];
	// halt("taper return pour continuer");
endfunction

