# Distributions de Kolmogorov (maple)
# 7 juillet 2010

with(LinearAlgebra):

# puissance de matrices avec conservation d'exposants

accuracy := 1e-8;

mPower := proc(A, eA, m, n)
	local B, eB;# double, int
	B := Matrix(m);
	if n = 1 then return A, eA;
	else
		(B, eB) := mPower(A, eA, m, floor(n/2));
		B := B.B; eB := 2*eB;
		if modp(n, 2) = 0 then return B, eB;
		else return A.B, eA;
		end if;
	end if;
end proc;

kolmogorovcdf := proc(d, n)
	local k, m, eH, eQ, i, j, g,# int
		p, h, s, H, Q;# double
	# Supprimer les prochaines lignes pour une
	# pr\'ecision sup\'erieure \`a 7 d\'ecimales
	s := d*d*n;
	if (s > 7.24 or (s > 3.76 and n > 99)) then
		return evalf(1-2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s));
	end if;
	k := floor((n*d)+1); m := 2*k-1; h := k-n*d;
	H := Matrix(m); Q := Matrix(m);
	for i from 1 to m do
		for j from 1 to m do
			if (i+1 < j) then H[i,j] := 0;
			else H[i,j] := 1;
			end if;
		end do;
	end do;
	for i from 1 to m do
	# le calcul des puissances $h^n$ pourrait \^etre am\'elior\'e
	# au niveau du temps. La question se pose au niveau
	# de la pr\'ecision...
		H[i,1] := H[i,1]-h^i;
		H[m,i] := H[m, i]-h^(m-i+1);
	end do;
	if 2*h-1 > 0 then H[m,1] := H[m,1]+(2*h-1)^m; end if;
	for i from 1 to m do
		for j from 1 to m do
			if i-j+1 > 0 then
				for g from 1 to (i-j+1) do
					H[i,j] := H[i,j]/g;
				end do;
			end if;
		end do;
	end do;
	eH := 0; (Q, eQ) := mPower(H, eH, m, n); p := Q[k,k];
	for i from 1 to n do
		p := p*i/n;
		if p < 1e-140 then
			p := p*1e140; eQ := eQ-140;
		end if;
	end do;
	return p*(10^eQ);
end proc;

kolmogorovicdf := proc(p, n)
	global accuracy;
	local a, b;# double
	if p > 1-accuracy then return 1;
	elif p < accuracy then return 0;
	else a := 0; b := 1;
		while b-a >= accuracy do
			if kolmogorovcdf((a+b)/2, n) < p then
			a := (a+b)/2;
			else b := (a+b)/2;
			end if;
		end do;
	return evalf((a+b)/2);
	end if;
end proc;

# pas encore teste, implementation asymptotique (Dudley)

klmcdf := proc(x, n)
	global accuracy;
	local z, p, term, k, pm;
	if x <= 0 then return 0;
	elif x >= 1 then return 1;
	else
		z := evalf(-2*(sqrt(n)+0.12+0.11/sqrt(n))^2*x*x);
		p := 1.0; k := 1; pm := -1; term := 1;
		while abs(term) > accuracy and k < 100 do
			term := evalf(2*pm*exp(z*k*k));
			p := evalf(p+term); pm := -pm; k := k+1;
		end do;
		return evalf(max(min(p, 1), 0));
	end if;
end proc;

klmicdf := proc(p, n)
	global accuracy;
	local a, b;
	if p > 1-accuracy then return 1;
	elif p < accuracy then return 0;
	else
		a := 0; b := 1;
		while b-a >= accuracy do
			if klmcdf((a+b)/2, n) < p then
				a := (a+b)/2;
			else b := (a+b)/2;
			end if;
		end do;
		return evalf((a+b)/2);
	end if;
end proc;

# 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 : F := x -> x; k := KSstat(u, F); p := 1-kolmogorovcdf(k, n);

KSstat := proc(u, F)
	local k, v, n, i;
	n := max(op(2, op(2,eval(u))), nops(u));
	v := sort(u);
	k := 0;
	for i from 1 to n do
	k := evalf(max(k, max(F(v[i])-(i-1)/n, i/n-F(v[i]))));
	end do;
	return k;
end proc;