(* Instructions and examples *)
(* 1. The use of this module requires the Numerical Calculus package: Needs["NumericalCalculus`"] *)
(* 2. Then it is necessary to define the following parameters *)
(* i)  Samples - is a vector with the samples size considered for each population *)
(* ii) type - can assume two values 1 or 2 *)
(*     if type=1 the output will be of the probability density function in a given value x *)(*     if type=2 the output will be of the cumulative distribution function in a given value x *)
(* iii) p - the number of variables in each of the populations *)
(* iv)  m - the number of exact moments matched by the approximating distribution *)
(* v)   x - the running value *)


Needs["NumericalCalculus`"];
EqMatrices[Samples_,type_,p_,m_,x_]:=Module[{(*q,NN,n,moments1,moments2,b,p1,r1,r2,l1,t,pp,ppp,pe,r1s*)},If[AllTrue[Samples,#>p&]&&IntegerQ[p]&&IntegerQ[m]&&AllTrue[Samples,IntegerQ]&&AnyTrue[{1,2},#==type&],ClearSystemCache;q=Length[Samples];NN=Total[Samples]-q;n=Samples-1;
xkl[k_,l_]:=Sum[n[[i]],{i,1,k-1}]+l;Clear[t];
moments1=SetPrecision[Table[(-1)^(-h)*ND[NN^(NN*t*p/2)/Product[n[[k]]^(n[[k]]*t*p/2),{k,1,q}]*Product[Gamma[(NN+1-j)/2]/Gamma[(NN+1-j)/2+NN/2*t]*Product[Gamma[(n[[k]]+1-j)/2+n[[k]]/2*t]/Gamma[(n[[k]]+1-j)/2],{k,1,q}],{j,1,p}],{t,h},0,Terms->7,WorkingPrecision->250,Scale->0.00000001],{h,1,If[m<4,4,m]}],750];
b=Flatten[Table[Table[Table[(NN+1-j)/(2*NN)+(xkl[k,l]-1)/NN-((n[[k]]+1-j)/(2*n[[k]])+(l-1)/n[[k]]),{l,1,n[[k]]}],{k,1,q}],{j,1,p}]];Clear[p1,r1,r2,l1];
{p1,r1,r2,l1}=Sort[Cases[{p1,r1,r2,l1}/.NSolve[Table[moments1[[h]]==Moment[MixtureDistribution[{p1,1-p1},{GammaDistribution[r1,1/l1],GammaDistribution[r2,1/l1]}],h],{h,1,4}],{p1,r1,r2,l1}],{_Real,_Real,_Real,_Real}]][[2]];
If[Abs[r1-Total[b]]<1,r1s=Total[b],r1s=r1];pe=Table[Unique["pe"],{m}];
moments2=SetPrecision[Table[Moment[MixtureDistribution[Flatten[{pe,{1-Total[pe]}}],Table[GammaDistribution[r1s+i,1/l1],{i,0,m}]],h],{h,1,m}],750];
ppp=Rationalize[Flatten[pe/.NSolve[Table[moments2[[h]]==moments1[[h]],{h,1,m}],pe]],0];
pp=Sum[ppp[[j]],{j,1,Length[ppp]}];l1=Rationalize[l1,0];r1s=Rationalize[r1s,0];
If[type==1,SetPrecision[Sum[ppp[[j+1]]*PDF[GammaDistribution[r1s+j,1/l1],x],{j,0,m-1}]+(1-pp)*PDF[GammaDistribution[r1s+m,1/l1],x],16],If[type==2,SetPrecision[Sum[ppp[[j+1]]*CDF[GammaDistribution[r1s+j,1/l1],x],{j,0,m-1}]+(1-pp)*CDF[GammaDistribution[r1s+m,1/l1],x],16]]],Print["Incorrect choice of parameters values"]]]

(* Example for the approximation of the cdf of negative logarithm of the likelihood ratio test statistic *)

Timing[SetPrecision[EqMatrices[{20,30,10,20},2,3,6,3/4],12]]

(* Example for the approximation of the pdf of negative logarithm of the likelihood ratio test statistic *)

Timing[SetPrecision[EqMatrices[{20,30,10,20},1,3,6,3/4],12]]

(* Example for the approximation of the cdf of likelihood ratio test statistic *)

Timing[SetPrecision[1-EqMatrices[{20,30,10,20},2,3,6,-Log[3/4]],12]]

(* Example for the approximation of the pdf of likelihood ratio test statistic *)

Timing[SetPrecision[EqMatrices[{20,30,10,20},1,3,6,-Log[3/4]]*4/3,12]]