This is a Maple program to compute the "dependant case" (PNFC2)
To use these programs you must use Maple software program. When Maple is executed, open the file
"PNFC2.mws".
Before using a function of the program, you have to compile all the procedure: put the
pointer on the first line of the first function (MakeSeq) and compile with the "enter" key. Continue with
"enter" until "end" at the end of program.
Continue with the same method to test examples
Sum of terms of a list L, from index n1 to index n2 (indexes begin from 1...)
>
SOMME := proc(L::list, n1::integer , n2:: integer) local i , longueur , total::integer;
longueur := nops(L) ; total:= 0 ;
if n2 > longueur or n1 > n2 or n1 <=0 then
total := 0 ;
else for i from n1 to n2 do
total := total + L[i] ;
od;
fi;
total;
end:
return the maximum index j of the list L such that the sum of the last terms is <= x -L[j]
TestMax := proc( L::list, k::integer , x::float)
local n , Splus , i, j1 ;
j1 := 1 ; Splus:= 0; n := nops(L) ;
Splus:= SOMME(L, n-k + 2, n) ;
while Splus <= x - L[j1] and j1 < n do j1 := j1 +1 ; od;
j1 := j1-1 ;
end:
return the minimum index j of the list L such that the sum of k-2 terms is <= x -L[j]
TestMin := proc (L::list , k::integer, x::float )
local Smoins , i, j2 ,n ;
j2 := 1 ; n := nops(L); Smoins:= 0 ;
Smoins := SOMME(L, j2 + 1 , j2 +k -1);
while j2 <=n and j2 + k -1 <=n -1 and Smoins <= x - L[j2] do
j2 := j2 + 1 ;
Smoins := Smoins - L[j2] + L[ j2 + k -1] ;
od;
j2 ;
end:
CasParticulier := proc (L::list, x:: float) local card:: integer , i :: integer , longueur :: integer;
i:= 1 ; card := 0; longueur:= nops(L);
while i <= longueur and L[i] <= x do
card := i ; i := i+1 ;
od;
card;
end:
create from the list L the new list of terms, the indexes of which are >= k.
TruncList := proc( L::list , k :: integer)
local L1 :: list , j :: integer , m :: integer;
m := nops(L);
if k <= m then
L1 := [seq(x[i],i=1.. m-k +1)];
for j from 1 to m-k+1 do
L1[j] := L[k+j-1] ;
od;
L1;
else print("Erreur de données entrées");
fi;
end:
>
Computation of max and min of sum of terms of lists.
>
Max1_function := proc(L1 ::list, ld :: float , k ::integer , l ::integer)
local q :: integer , M :: float , L2 :: list ;
q := nops(L1); M:= 0 ; L2 := ld * L1 ;
if k = 1 then
M := SOMME (L1, q-l+1, q);
else
M := M + SOMME( L1, q-k+2, q) + SOMME (L2, q-k-l+2, q-k+1) ;
fi;
end:
>
MaxLamda_function := proc(L1 ::list, ld :: float , k ::integer , l ::integer)
local q :: integer , M :: float , L2 :: list ;
q := nops(L1); M:= 0 ; L2 := ld * L1 ;
if k = 1 then
M := SOMME (L1, q-k+1, q);
else
M := M + SOMME( L1, q-k+1, q) + SOMME (L2, q-k-l+2, q-k+1) ;
fi;
end:
>
Min_function := proc(L1 ::list, ld:: float, k ::integer , l ::integer, j :: integer)
local q :: integer , M :: float , L2:: list ;
q := nops(L1) ; M:= 0 ; L2 := ld*L1 ;
if j <= q - k -l + 1 then
M := M + SOMME (L1, j+1, j+k-1) + SOMME (L2, j+k, j+k+l-1);
fi;
M;
end:
Partial computations for favourable cases...
C_function := proc(M :: integer, n:: integer, k ::integer , l ::integer)
option remember ;
local j :: integer , C :: integer ;
C := 0 ;
for j from 1 to M do
C := C + binomial (n-j, k) * binomial (n-j-k, l-1 );
od;
end:
Computation for "lambda" = 1.
Calcul_NFC1 := proc ( L:: list, k :: integer , x ::float )
option remember ;
local i :: integer, j1 :: integer, j2:: integer , Nk :: integer , n :: integer , Lo :: list ;
Lo:= L; n := nops(Lo) ; Nk := 0 ;
if k = 0 then Nk:= Nk +1 ; else
if k > n then Nk else
if k = 1 then Nk := Nk + CasParticulier (Lo, x); else
j1 := TestMax( Lo, k , x);
j2 := TestMin(Lo, k , x );
if j1 > 0 then
for i from 1 to j1 do
Nk := Nk + binomial (n-i, k-1) ;
Nk;
od;
fi;
if j1 + 1 <= j2 -1 then
for i from j1 +1 to j2 -1 do
Nk := Nk + Calcul_NFC1 ( TruncList(Lo, i+1 , n), k-1 , x - Lo[i] ) ;
od;
fi;
fi;
fi;
fi;
if Nk = (NULL) then Nk := 0 fi;
Nk ;
end:
>
The main computation, independant case... function NFC
L: the sequence ; "ld"= lambda, a real number between 0 and 1 ;
k,l = integers such that k+l <= n ; x a real number.
>
NFC := proc ( L:: list, k :: integer , l ::integer, ld :: float , x :: float)
option remember
;local ldL :: list , N :: integer , MaxLd :: integer , n :: integer, m1 :: integer ,
j:: integer , M1 :: integer , Mld :: integer , mld :: integer, Max1 , min1, minLd, S;
ldL := ld*L ; n := nops( L);
n= size of the list ;
ld = coeff "lambda" ;
x = a real number ; k,l = integers such that k+l <= n ;
the function returns : N = Nk = number of favourable cases ... <= x
N:= 0 ;
if x >= 0 then
if k+l <= n then
if k = 0 and l = 0 then N := 1
else
if k = 0 then
N := N + Calcul_NFC1 ( ldL, l, x);
fi;
if l = 0 then
N := N + Calcul_NFC1 ( L, k, x);
fi;
fi;
Begining
>
if k >0 and l >0 then
détermination of Mld
MaxLd := MaxLamda_function (L, ld, k, l);
Mld := 1 ;
while MaxLd <= n-k-l and MaxLd <= x - ld * L[Mld] do
Mld := Mld+1 ;
od;
Mld := Mld - 1 ;
détermination of M1
Max1 := Max1_function ( L , ld , k, l);
M1 := Mld -1 ;
while M1 <= n-k-l and M1 >0 and Max1 <= x - L[M1] do
M1 := M1+1;
od;
if M1 >1 then M1 := M1-1 else M1 := 0; fi;
détermination of m1
m1 := M1 +1 ;
if m1 <= n -k -l +1 and m1 < n then
min1 := Min_function ( L, ld, k , l , m1 ) ;
while m1 <= n - k -l and min1 <= x - L[m1] do
min1 := min1 + L[m1 + k] -L[m1+1] + ld * ( L[m1+k+l] - L[m1 + k]) ;
m1 := m1+1 ;
od;
else m1 := n - k - l + 2 ;
fi;
détermination of mld = m_lambda
mld := m1;
if mld <= n -k -l +1 then
minLd := Min_function ( L, ld, k+1 , l , mld ) ;
while mld+ k +2 <= n and minLd <= x - ld * L[m1] and mld <= n - k -l and mld >0 do
minLd := minLd + L[mld + k + 1] -L[ mld+1] + ld * ( L[mld+k+2] - L[mld + k + l]) ;
mld := mld+1 ;
od;
else mld := n - k - l + 2 ;
fi;
mld := mld + 1 ;
m1 := m1 + 1 ;
The computation...
N := N + C_function( M1 , n , k-1, l +1 ) + C_function( Mld , n , k , l ) ;
if M1 + 1 <= m1 -1 then
for j from M1 + 1 to m1 -1 do
if n >= 1 then S := TruncList(L, j+1);
N := N + NFC( S , k-1 , l, ld , x - L[j]) ;
fi
od;
fi;
if Mld+1 <= mld -1 then
for j from Mld +1 to mld -1 do
if n >= 1 then S := TruncList(L, j+1);
N := N + NFC( S , k , l-1, ld , x - ld * L[j]) ;
fi;
od;
fi;
fi;
fi;
else ;;
fi;
N; N - (NULL);
end:
Make a sequence of "n" random positive integers between 0 and "max"
>
MakeSeq := proc (n :: integer, max :: integer)
local j :: integer, S, L ;
S := seq(rand(0..max)(),i=1..n):
L := [seq(x[i],i=1.. n)];
for j from 1 to n do
L[j] := S[j] :
od:
L:=sort(L);
end:
>
>
The function P_NFC gives the probabilities for dependant case: lambda = lambda1 / lambda2
S a sequence ; u1 = u_0 (i) ; u2 = u_1 (j).
>
>
PNFC2 := proc ( S :: list, lambda1 :: integer, lambda2 :: integer , u1 , u2 , t :: float)
local NTk:: float, k :: integer , n ,mu_k, M1 :: integer, m1 :: integer ;
n := nops(S) ;
M1 := min( u1, u2) : m1 := max ( 0, u1 + u2 - n );
NTk := 0 ;
for k from m1 to M1 do
mu_k := binomial (u1, k)* binomial(n- u1, n-u1 -u2 +k)/( binomial(n,u2) * binomial (n,k) *binomial(n-k, n - u1 - u2 +k)) ;
NTk := NTk + mu_k * NFC(S, k, n - u1 - u2 +k , 0.5 , t) ;
od;
NTk; convert( NTk - (NULL), float) ;
end:
An example:
S := [1 , 0.6 , 0.8 1 , 0.34 , 0.26 , 0.448 , 0.484 , 0.18 , 0.364 , 0.248 , 0.196 , 0.212 , 0.228]
> S := [1 , 0.6 , 0.8 , 1 , 0.34 , 0.26 , 0.448 , 0.484 , 0.18 , 0.364 , 0.248 , 0.196 , 0.212 , 0.228];
Sort the list and make it an integer list by multiplying by 1000.
> S:= sort(S); S:= 1000*S;
> S0 := [180, 196, 212, 228, 248, 260, 340, 364, 448, 484, 600, 800, 1000, 1000];
X= 3 (x1000) (with a dot after 3000 to force a real number); lambda = 1/2, u1=5 ; u2=7;
> tm := time():PNFC2 ( S0, 1 , 2 , 5, 7 , 3000.); time() -tm;
> tm := time():PNFC2 ( S0, 1 , 2 , 5, 10 , 3000.);time() -tm;
>