# restart; read `c:/Users/Asus/Google Drive/Aek/
# MomentRamsey/Ramsey.txt`;
# read `c:/Users/Asus/Google Drive/Guess.txt`;

with(combinat)[stirling1]:

# March 30, 2018
##################################
# Section 1: 
# Moment of Delaporte distribution
# E[X^k] and E[(X-mu)^k]
##################################

# Moment of Delaporte distribution
# Output: E[X^m]
# Try: seq(MoDe(m,a,b,s),m=1..4);
MoDe := proc(m,a,b,s) option remember;
local t,A;

if m = 0 then 
	return(1);
fi:	

A :=exp(s*(exp(t)-1))/(1-b*(exp(t)-1))^a;
expand(subs(t=0,diff(A,t$m)));
end:


# Moment about the mean of Delaporte distribution
# Output: E[(X-mu)^m]
# Try: seq(MoMeanDe(m,a,b,s),m=1..6);
# sort(MoMeanDe(10,a,b,0));
MoMeanDe := proc(m,a,b,s) option remember;
local i,F;

F := add( (-1)^(i)*binomial(m,i)*
	MoDe(m-i,a,b,s)*MoDe(1,a,b,s)^i ,i=0..m);

return(expand(F));
end:


# Falling Moment of Delaporte dist.
# using probability generating function
# Output: E[(X)_m]
# # Try: FallMoDe(1,a,b,s);
FallMoDe := proc(m,a,b,s) option remember;
local z,A;

if m = 0 then 
	return(1);
fi:	

A :=exp(s*(z-1))/(1-b*(z-1))^a;
expand(subs(z=1,diff(A,z$m)));
end:

# Falling Moment of Delaporte dist.
# Output: E[(X)_m]
# Try: FallMoDe2(1,a,b,s);
# sort(FallMoDe2(6,a,b,s),s);
FallMoDe2 := proc(m,a,b,s) option remember;

add(stirling1(m,i)*MoDe(i,a,b,s),i=0..m);

end:



##################################
# Section 2: 
# Coeff of leading term n^L of E[(X-mu)^t]
##################################


# Section 2.1: E[(X-mu)^2] with 
# computation and general computation


# Check coeff of n^(2*k-1) of E[(X-mu)^2]
# = E[X^2]-E[X]^2.
# Try: (1/(2^bino(k,2)-1)/k!)^2*First2(k);
First2 := proc(k) option remember;
local A;

A := Stir1(2*k,1)+k^2*Stir1(2*k-1,0)
-add(Stir1(k,i)*Stir1(k,1-i),i=0..1);

factor(A);
end:

# Check coeff of n^(2*k-2) of E[(X-mu)^2]
# Try: (1/(2^bino(k,2)-1)/k!)^2*Second2(k);
Second2 := proc(k) option remember;
local A;

A := Stir1(2*k,2)+k^2*Stir1(2*k-1,1)
+k^2*(k-1)^2/2*Stir1(2*k-2,0)
- add(Stir1(k,i)*Stir1(k,2-i),i=0..2);

factor(A);
end:

# Check coeff of n^(2*k-3) of E[(X-mu)^2]
# Try: (1/(2^bino(k,2)-1)/k!)^2*Third2(k);
Third2 := proc(k) option remember;
local A;

A := Stir1(2*k,3)+k^2*Stir1(2*k-1,2)
+1/2*(k*(k-1))^2*Stir1(2*k-2,1)
+4/3!*(k*(k-1)*(k-2))^2*Stir1(2*k-3,0)
- add(Stir1(k,i)*Stir1(k,3-i),i=0..3);

factor(A);
end:


# General formula for coeff of n^(2*k-t) of E[(X-mu)^2]
# Try: seq(GenTerm2(t,k),t=1..4);
# First2(k); Second2(k); Third2(k);
GenTerm2 := proc(t,k) option remember;
local i,C,A;

#C:= 2/2^binom(k,2)/k!:

A:= add(Factor2(i,k)*Stir1(2*k-i,t-i),i=0..t)
- add(Stir1(k,i)*Stir1(k,t-i),i=0..t);

factor(C^2*A);
end:

# Factor2= Term/C, where term is more like
# a probability of mono / k!/s!
# Try: seq(Factor2(s,k),s=0..5);
Factor2 := proc(s,k) option remember;
local con;

if s=0 or s=1 then
	con:=1;
else
	con :=2^binomial(s,2)/2; 
fi:	

con*Poch(k,s)^2/s!;

end:

######################################


# Section 2.2: E[(X-mu)^3] with 
# computation and general computation


# Check coeff of n^(3*k-1) of E[(X-mu)^3]
# = E[X^3]-3*E[X^2]*E[X]+2*E[X]^3

# Try: (1/(2^bino(k,2)-1)/k!)^3*First3(k);
First3 := proc(k) option remember;
local i,j,A;

A := Stir1(3*k,1)+3*k^2*Stir1(3*k-1,0)
-3*(Stir1(2*k,0)*Stir1(k,1)
+(Stir1(2*k,1)+k^2*Stir1(2*k-1,0))*Stir1(k,0))
+2*add(add(Stir1(k,i)*Stir1(k,j)*Stir1(k,1-i-j)
,j=0..1-i),i=0..1);

factor(A);
end:


# Check coeff of n^(3*k-2) of E[(X-mu)^3]
# Try: (1/(2^bino(k,2)-1)/k!)^3*Second3(k);
Second3 := proc(k) option remember;
local i,j,A;

A := Stir1(3*k,2)+3*k^2*Stir1(3*k-1,1)
+k^3*Stir1(3*k-2,0)+3/2*(k*(k-1))^2*Stir1(3*k-2,0)
+3*k^3*(k-1)*Stir1(3*k-2,0)
-3*(Stir1(2*k,0)*Stir1(k,2)
+(Stir1(2*k,1)+k^2*Stir1(2*k-1,0))*Stir1(k,1)
+(Stir1(2*k,2)+k^2*Stir1(2*k-1,1)
+1/2*(k*(k-1))^2*Stir1(2*k-2,0))*Stir1(k,0))
+2*add(add(Stir1(k,i)*Stir1(k,j)*Stir1(k,2-i-j)
,j=0..2-i),i=0..2);

factor(A);
end:


# General formula for coeff of n^(3*k-t) 
# of E[(X-mu)^3]

# Try: seq(GenTerm3(t,k),t=1..4);
# First3(k); Second3(k);
GenTerm3 := proc(t,k) option remember;
local i,j,s,f,C,A;

#C:= 2/2^binom(k,2)/k!:

A := add(add(Factor3(f,k),f in Config3(i))
*Stir1(3*k-i,t-i),i=0..t)
-3*add(add(MyFactor(i,k)*Stir1(2*k-i,s-i),i=0..s)
*Stir1(k,t-s),s=0..t)
+2*add(add(Stir1(k,i)*Stir1(k,j)*Stir1(k,t-i-j)
,j=0..t-i),i=0..t);


factor(C^3*A);
end:

# Output: set of [a,b,c,d] where 2a+b+c+d=t
# Try: Config3(3);
Config3 := proc(t) option remember;
local a,b,c,T;

T := {};

for a from 0 to t/2 do
	for b from 0 to t-2*a do
		for c from 0 to t-2*a-b do
			T := T union {[a,b,c,t-2*a-b-c]};
od:od:od:

T;
end:


# Factor of each term (Fac= Term/C)
# Try: Factor3([1,0,0,0],k);
Factor3 := proc(CF,k) option remember;
local i,S,T,score,con;

S := [CF[1]+CF[2]+CF[3], 
CF[1]+CF[2]+CF[4], CF[1]+CF[3]+CF[4]];

T := [CF[1]+CF[2], CF[1]+CF[3], CF[1]+CF[4]];
score := 0;

for i from 1 to 3 do
	if T[i] >=2 then
		score := score+1;
	fi:
od:	

if score >=2 then
	con := 1/4;
elif score =1 then
	con := 1/2;
else
	con := 1;
fi:	

for i from 1 to 3 do
		con :=con*2^binomial(T[i],2); 
od:

con := con/2^binomial(CF[1],2);	
	
con*Poch(k,S[1])*Poch(k,S[2])*Poch(k,S[3])
/CF[1]!/CF[2]!/CF[3]!/CF[4]!;

end:


######################################

# Section 2.3:
# Helper functions

# Stirling formula of the first kind
# by fitting the polynomial
# Try: seq(Stir1Pol(k,s),s=1..4); 
Stir1Pol := proc(k,s) option remember;
local i,A,x;

A := [seq(stirling1(i,i-s),i=s..3*s+3)];
subs(x=k,GuessPol(A,s,x));
		
end:

# binomial with symbolic n
# Try: binom(n,3);
binom := proc(n,k) option remember;
local i;

mul(n-i,i=0..k-1)/k!;
end:

# Pochammer
# Try: Poch(n,3);
Poch := proc(n,k) option remember;
local i;

mul(n-i,i=0..k-1);
end:


##################################
# Section 3: 
# Converting Bonfer to Straight Moment
##################################

# Try: ConBonfer(A,1);
ConBonfer := proc(A,m) option remember;
local i,s;

add((-1)^s*add(A[i]*Stir1(s,s-i),i=0..s)/s!,s=0..m);

end:

# Falling Moment of Poisson dist
# Try: [seq(FallMoPoi(s,m),m=0..7)];
FallMoPoi := proc(s,m) option remember;
local z;

simplify(subs(z=1,diff(exp((z-1)*s),z$m)));
end:



##############################
# Section 4: 
# Verify Delaporte for small k
# in R(k,k)
##############################


##################################
# Section 4.1: 
# Use my moment (with least squares
# method) to approximate alpha, 
# beta, lambda
##################################

# I try to fit a,b,lambda by the
# fitting the least squares of
# the moments. But it did not give
# a good fit anyway. (compare to
# Aaron's result on page 11)
# Try: FitD(4,14);
FitD := proc(k,n) option remember;
local i,K,N,Mea,MM2,MM3,MM4,MM5
,A,a,b,s,eq,var;

Mea := subs({N=n},Moment(1,k,N,2));
MM2 := subs({N=n},SecondMomentMean(k,N,2));
MM3 := subs({N=n},ThirdMomentMean(k,N,2));
MM4 := subs({N=n},FourthMomentMean(k,N,2));
MM5 := subs({N=n},FifthMomentMean(k,N,2));

A := [MoDe(1,a,b,s)-Mea,MoMeanDe(2,a,b,s)-MM2,
MoMeanDe(3,a,b,s)-MM3, MoMeanDe(4,a,b,s)-MM4
, MoMeanDe(5,a,b,s)-MM5];

A := add(A[i]^2,i=1..nops(A));

eq := {diff(A,a),diff(A,b),diff(A,s)};
var := {a,b,s};

evalf(solve(eq,var));
end:



##################################
# Section 4.2: 
# Given alpha, beta, lambda, compare
# moment from delaporte and actual
##################################

# Try: Example 1 from Aaron's paper page 11. 
# TestAa(4,14,1.84,7.89,16.75);
# Second example
# TestAa(4,20,3.74,15.46,93.57);
# Third example
# TestAa(5,49,9.45,163.51,2178.57);

TestAa := proc(k,n,a,b,s) option remember;

print("Delaporte distribution");
print("Mean", MoDe(1,a,b,s));
print("Variance", MoMeanDe(2,a,b,s));
print("Third Moment about mean", MoMeanDe(3,a,b,s));
print("Fourth Moment about mean", MoMeanDe(4,a,b,s));
print("Fifth Moment about mean", MoMeanDe(5,a,b,s));
print("");
print("Theorectical value");
print("Mean", evalf(Moment(1,k,n,2)));
print("Variance", evalf(SecondMomentMean(k,n,2)));
print("Third Moment about mean", 
evalf(ThirdMomentMean(k,n,2)));
print("Fourth Moment about mean", 
evalf(FourthMomentMean(k,n,2)));
print("Fifth Moment about mean", 
evalf(FifthMomentMean(k,n,2)));
end:

# How about checking Bonferoni from 
# the numbers a,b,lambda?
# Does not converge yet. need very
# big moment
# Output: Check 1-E[X], 
# 1-E[X]+E[(X,2)]-E[(X,3)], ...
# Try: 
# Example 1
# CheckBon(1.84,7.89,16.75,1000);
# Second example
# CheckBon(3.74,15.46,93.57,40);
# Third example
# CheckBon(9.45,163.51,2178.57,20);
CheckBon := proc(a,b,s,L) option remember;
local i,m,T;

T := 0;

for i from 0 to L do
	T := T+add((-1)^m*FallMoDe(m,a,b,s)/m!,m=2*i..2*i+1);
	print("Term", 2*i+1, "Value", T);
od:

return();
end: