###############################################
# Filename: /stat/chien/JULY95/MOMENT/maplecode
# Date:     July 31, 1996
# By:       Fred W. Huffer, Florida State Univ.
#           Chien-Tai Lin, Tamkang Univ.
# Codes for Tables in Huffer and Lin (1995).  
#
#    read `/stat/chien/JULY95/MOMENT/maplecode`;
#
#  Type this READ command to read codes for all of 
#  the methods. You may change the directory to 
#  your own directory which contains this file, eg., 
#         read `directory/maplecode`;
#
#  To execute the readmom procedure next, you may
#  also want to change the directory inside the
#  READ command to your own directoy that contains
#  the files which have the polynomial outputs for
#  different values of m, eg.,
#        read `directory/filename-`.mval;  
#
# Tables 1 and 3:
#    readmom(mval); (Enter the values of n,m,d.)
#    setval(dval,nval);
# (Variables n,m,d,w=n-m+1 are global.)
#
#    Table 1:
#      mc2();
#      lp4(0); ( or lp4(1); to print out
#               the values of lambda's.)   
#      cpg4(1); (jmax=1)
#      klb();
#      bounds();
#
#    Table 3;
#      cpg2(5); (jmax=5)
#      cpg4(5);
#      gcp(5); 
#
# Table2:
#    fsmval(mval,dval,nval);
#    cpg2(1);  
#    mc2();
###############################################

Digits:=20;

readmom:=proc(mval)
######################################################################
# This procedure reads in the expressions for the first four moments
#    which are named m1, m2, m3, m4.
# Example:  To get the expressions for m=7, execute readmom(7).
######################################################################
global m,mmmmm,m1,m2,m3,m4;
m:=mval;
mmmmm:=mval; 
read `/stat/chien/JULY95/MOMENT/work1-`.mval;
#################################################################
# Read the output files from your own directory inside the codes.
#################################################################
end;

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

setval:=proc(dval,nval)
##########################################################################
# This computes the first four moments and first four cumulants 
#    for d=dval and n=nval.
# These are named mom1, mom2, mom3, mom4, c1, c2, c3, c4.
# You must execute the procedure readmom before using setval
#    in order to specify the value of m and read in the expressions
#    for the moments.  However, as long as you stick with the same
#    value of m, you only need to execute readmom once.  You can then run 
#    setval as many times as you like.  
# You may wish to set Digits before using setval.
##########################################################################
global m,mmmmm,d,ddddd,n,nnnnn,w,wwwww,
       c1,c2,c3,c4, mom1,mom2,mom3,mom4;
d:=dval;
n:=nval;
ddddd:=dval;
nnnnn:=nval;
wwwww:=nnnnn-mmmmm+1;
w:=wwwww;
mom1:=evalf(m1);
mom2:=evalf(m2);
mom3:=evalf(m3);
mom4:=evalf(m4);
c1:=mom1;
c2:=mom2-mom1^2;
c3:=mom3-2*mom1*c2-mom2*mom1;
c4:=mom4-3*mom1*c3-3*mom2*c2-mom3*mom1;
NULL;
end;

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

b:=proc(j,k)
  Heaviside(n-j-k)*binomial(n-j,k);
end;

R:=proc(s,t)
   if 1-t*d<=0 then 
      0;
   else
      binomial(n,s)*d^s*(1-t*d)^(n-s);
   fi;
end;

##############################################################
# Evaluate the first two moments from G and F for larger m.
# Faster Routine to get the first two moments
##############################################################

acc2:=proc(a,dim)
local i,j,tot;
for j from 1 to dim do
tot:=0;
for i from 0 to j-1 do
tot:=tot+a[i,j];
a[i,j]:=a[i,j-1]+tot;
od;
a[j,j]:=a[j-1,j-1]+2*tot+a[j,j]
od;
eval(a);
end;


acc1:=proc(a,dim)
local j;
for j from 1 to dim do
a[j]:=a[j-1]+a[j];
od;
eval(a);
end;

####################
# The first moment 
####################

onemom:=proc(mm,dd,nn)
local y,i,g,ans;
y:=1-dd;
if y<=0 then RETURN(`Failed.`); fi;
g:=array(0..(mm-2),[seq(evalf(binomial(nn,i)*dd^i*y^(nn-i)),i=0..(mm-2))]);
acc1(g,mm-2);
ans:=(nn-mm+1)*(1-g[mm-2]);
ans;
end;

####################
# The second moment
####################

twomom:=proc(mm,dd,nn)
local f0,i,j,y,z,g,f,ans;
y:=1-dd;
z:=1-2*dd;
if y <= 0 or z <=0 then RETURN(`Failed.`); fi;
g:=array(0..(mm-2),[seq(evalf(binomial(nn,i)*dd^i*y^(nn-i)),i=0..(mm-2))]);
f:=array(0..(mm-2),0..(mm-2));
for j from 0 to (mm-2) do for i from 0 to j do
f[i,j] := evalf(binomial(nn,i)*binomial(nn-i,j)*d^(i+j)*z^(nn-i-j));
od;
od;
ans:=0;
acc1(g,mm-2);
ans:=(nn-mm+1)*(1-g[mm-2])+(nn-mm+1)*(nn-mm)*(1-2*g[mm-2]);
acc1(g,mm-3);
acc1(g,mm-3);
ans:=ans+4*(nn-mm)*g[mm-3];
acc1(g,mm-4);
acc1(g,mm-4);
ans:=ans-12*g[mm-4];
acc2(f,mm-2);
ans:=ans+(nn-2*mm+3)*(nn-2*mm+2)*f[mm-2,mm-2];
acc2(f,mm-3);
ans:=ans-2*(nn-2*mm+3)*f[mm-3,mm-3];
acc2(f,mm-4);
ans:=ans+2*f[mm-4,mm-4];
ans;
end;

###########################################################
# This procedure computes the first two moments and thus
# the first two cumulants by faster routine in order to 
# get the result for larger m=mval,d=dval,n=nval.
# Run mc2() and cpg2(1) procedures after this procedure can 
# obtain the values of Table 2.
############################################################
fsmval:=proc(mval,dval,nval)
   global n,m,d,w,mom1,mom2,c1,c2;
   n:=nval;
   d:=dval;
   m:=mval;
   w:=n-m+1;
   if m<=3 then ERROR(`Try bigger m`); fi;
   mom1:=onemom(m,d,n);
   if type(mom1,string) then ERROR(`Try another d.`); fi;
   if n>=2*(m-1) then
      mom2:=twomom(m,d,n);
      if type(mom2,string) then ERROR(`Try another d.`); fi;
   else
      ERROR(`Failed. n<2(m-1). Try another n or m.`);
   fi;
   c1:=mom1;
   c2:=mom2-mom1^2;
   NULL;
end;

max_val_i:=30;

######################################################
# BOUNDS: max_val_i can be changed if you like.
######################################################
bounds:=proc()
local i,j,p,ps,lb,ub,y;
i:='i';
j:='j';
p:=collect(simplify(expand(1-(y-i)*(y-i-1)*(y-j)*(y-j-1)/(i*(i+1)*j*(j+1)))),y);
ps:=unapply(subs(y^4=mom4,y^3=mom3,y^2=mom2,y=mom1,p),i,j);
lb:=max(seq(seq(ps(i,j),j=(i+2)..max_val_i+2),i=1..max_val_i));
i:='i';
p:=collect(simplify(expand(1-(y-1)*(y-w)*(y-i)*(y-i-1)/(w*i*(i+1)))),y);
ps:=unapply(subs(y^4=mom4,y^3=mom3,y^2=mom2,y=mom1,p),i);
ub:=min(seq(ps(i),i=2..min(max_val_i,w-2)));
printf(`LB=%.9f\n`,lb);
printf(`UB=%.9f\n`,ub);
end;

###########################################
#   MC2 : 0 < pi/s < 1 and (w-1)^2-4*c > 0 
###########################################
mc2 := proc()
local MC2,pi,c,s;
pi := mom1/w ;
c :=(c2 - w*pi*(1-pi))/(2*pi*(1-pi));
s := (1/2)*(w + 1 - ( (w-1)^2 - 4*c)^(1/2));
MC2 := 1 - (1-pi)*(1-pi/s)^(w-1);
if pi/s > 1 or pi/s <=0 then ERROR(`Failed on 0 < pi/s < 1.`); fi;
if (w-1)^2-4*c <=0 then ERROR(`Failed on (w-1)^2-4*c > 0`); fi;
printf(`MC2=%.9f\n`,MC2);
end;

###############################
#  CPG2 : c2-m1>0 and 0<=r<1
###############################
cpg2:=proc(jmax)
local r,i,j;
if c2-mom1 <0  then ERROR(`Failed on c2-m1 >= 0`); fi; 
r:=((c2/c1)-1)/((c2/c1)+1);
if r>1 or r< 0 then ERROR(`Failed on 0<= r < 1`); fi;
p.0 := exp(- 2*mom1/(1 + c2/mom1));
l.1:=c1*(1-r)^2;
for j from 2 to jmax do l.j:=l.1*r^(j-1); od;
for j from 1 to jmax do 
    p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,m))],`+`)/j;
od;
printf(`===CPG2===\n`);
pge.0:=1;
for j from 1 to jmax do 
    pge.j := pge.(j-1)-p.(j-1);
    printf(`%d,%.9f\n`,j,pge.j);
od;
end;

#################################################################
# CPG4 approximation
# all lambdas should be positive and 0<=x<1 (r in the paper) .
##################################################################
cpg4:= proc(jmax)
local f,v,F,a,xi,i,j,soln,g1,g2,g3,g4,r;
f:=x/(1-x);
with(linalg):
v:=vector(4);
v[1]:=simplify(x*diff(f,x));
for j from 2 to 4 do v[j] :=simplify(x*diff(v[j-1],x)) od;
F:=matrix(4,3,[1,2,v[1],1,4,v[2],1,8,v[3],1,16,v[4]]);
a:=matrix(3,1,[a0,a1,a2]);
xi:=multiply(F,a);
g1 := xi[1,1]=c1:
g2 := xi[2,1]=c2:
g3 := xi[3,1]=c3:
g4 := xi[4,1]=c4:
soln:= fsolve({g1,g2,g3,g4},{a0,a1,a2,x});
l.1:=evalf(subs(soln,a0+a2*x));
if l.1 < 0 then ERROR(`Failed on lambda 1.`); fi;
l.2:=evalf(subs(soln,a1+a2*x^2));
if l.2 < 0 then ERROR(`Failed on lambda 2.`); fi; 
l.3:=evalf(subs(soln,a2*x^3));
if l.3 < 0 then ERROR(`Failed on lambda 3.`); fi;
r:=evalf(subs(soln,x));
if r  < 0 or r > 1 then ERROR(`Failed on r.`); fi;
p.0:=exp(-evalf(subs(soln,a0+a1+a2*x/(1-x))));
for j from 4 to jmax do l.j:=l.3*r^(j-3); od;
for j from 1 to jmax do 
    p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,m))],`+`)/j;
od;
printf(`===CPG4===\n`);
pge.0:=1;
for j from 1 to jmax do 
    pge.j := pge.(j-1)-p.(j-1);
    printf(`%d,%.9f\n`,j,pge.j);
od;
end;

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

lpmin:=proc(q,deg,output)
##########################################################################
# This procedure is used by lp4.  It computes the minimum probability p
#   assuming a compound Poisson distribution for Y in which lambda[i]=L[i]
#   = zero for i>q and the nonzero lambda's satisfy constraints
#   of order deg.
#     deg=0 assumes only that L[i]>=0 for all i.
#     deg=1 assumes in addition that L[i]>=L[i+1] for all i.
#      (in other words, the first differences (L[i]-L[i+1]) are >= 0.)
#     deg=2 assumes in addition that the second differences are >=0.
##########################################################################
global c1,c2,c3,c4;
local kmom,i,j,e,constraint,func,L,A,T,tri,B,solmin,pmin;
kmom:=4;
A:=array(1..(kmom+1),1..q,[seq([seq(j^i,j=1..q)],i=0..kmom)]);
e:=array(1..q,1..1);
if deg=0 then
T:=&*();
else
# Create an upper triangular matrix of 1's (Is there a better way?)
tri:=array(1..q,1..q,[seq([seq((1+signum(j-i+.5))/2,j=1..q)],i=1..q)]);
T:=evalm(tri^deg);
fi;
B:=evalm(A &* T);
constraint:={seq(c.i=convert( [seq(B[i+1,j]*e[j,1],j=1..q)] ,`+`),i=1..kmom)};
func:=convert([seq(B[1,j]*e[j,1],j=1..q)],`+`);
solmin:=simplex[minimize](func,constraint,NONNEGATIVE);
if nops(solmin)=0 then
RETURN(`Failed.`); 
fi;
pmin:=1-evalf(subs(solmin,exp(-func))):
if output=1 then
assign(eval(solmin,1));
L:=evalm(T &* e);
print(`Lambda's are`,L);
fi;
pmin;
end;

###########
###   LP4
###########
lp4:=proc(output)
#######################################################################
# This procedure computes the LP4 approximation described in the paper.
# If output=1, the lambda's will be printed.
# It prints out the values of p, deg, and q.
# See further comments in the procedure lpmin.
#######################################################################
local deg,p,qmax,qqqqq,qend;
deg := 2;
qqqqq:=mmmmm+1;
qmax:=max(2*mmmmm,10);
while deg >= 0 do
   p := lpmin(qqqqq,deg,output);
   qend:=qqqqq;
   if type(p,string) and qend <=qmax  then 
      qqqqq:=qqqqq+5;
   elif qend>qmax then
      qqqqq:=mmmmm+1;
      deg:=deg-1;
   else 
      break 
   fi;
od;
if deg<0 then
printf(`No solution exists. Try a value of q > %d\n`,qend);
else
   printf(`lp4=%.9f,deg=%d,q=%d\n`,p,deg,qqqqq); 
fi;
end;

###########
#    GCP
###########
gcp := proc(jmax)
local pi,EY,mm,EI1I2,p,j,i;
mm:=m-1:
pi:=1-evalf(Sum('R(j,1)','j'=0..(mm-1)));
EY:=(n-mm)*pi;
EI1I2:=1-2*evalf(Sum('R(mm-2*j-1,1)','j'=0..floor((mm-1)/2)))
       +(-1)^(mm+1)*R(0,2);
p:=EI1I2/pi;
for j from 2 to mm do l.j := (n-mm)*pi*(1-p)^2*p^(j-1); od;
l.1:=EY-evalf(Sum('j*l.j','j'=2..mm));
p.0:=exp(-evalf(Sum('l.j','j'=1..mm)));
for j from 1 to jmax do 
    p.j:= convert([seq(i*(l.i)*(p.(j-i)),i=1..min(j,mm))],`+`)/j;
od;
printf(`===GCP===\n`);
pge.0:=1;
for j from 1 to jmax do 
    pge.j := pge.(j-1)-p.(j-1);
    printf(`%d,%.9f\n`,j,pge.j);
od;
end;

###############
#  KLB 
###############
klb:=proc()
local r,s2,lb;
s2:=(mom2-mom1)/2;
r:=floor(2*s2/mom1)+1;
lb:=2*(r*mom1-s2)/(r*(r+1));
printf(`klb=%.9f\n`,lb);
end;

###############
####   AP1-4
###############
ap:= proc()
  global c1,c2,c3,c4;
  local i,f,A,B,C,l;
  for i to 4 do
     f := (j,k) -> k^j:
     A := matrix(i,i,f);
     B := inverse(A);
     if i=1 then
        C := multiply(B,[c1]);
     elif i=2 then 
        C := multiply(B,[c1,c2]);
     elif i=3 then
        C := multiply(B,[c1,c2,c3]);
     else
        C := multiply(B,[c1,c2,c3,c4]);
     fi;        
     print(C);
     for l to i do  
         if C[l] < 0 then
              l := 6; 
         fi;
     od;
     if(l<=5) then 
        print(`AP`,i,1-exp(-sum('C[k]','k'=1..i)));
     else 
        print(`AP`,i);
     fi;             
  od;
end;