(*****************************************************************************************************************)
(* Instructions *)
(* 1. Please, first download the codes for the probability density (pdf) or the cumulative distribution (cdf) functions of a random variable with a GNIG distribtuion *)
(* 2. These codes are available in the web-page: https://sites.google.com/site/nearexactdistributions/ *)
(* 3. From the web-page you need to download the codes of:  *)
(*    a) GNIGpdf[r_,b_,l_,a_,w_] - for the pdf of the GNIG distribution *)
(*    b) GNIGcdf[r_,b_,l_,a_,w_] - for the cdf of the GNIG distribution *)
(*    c) Makec[r_,l_,p_] - auxiliary module *)
(* 4. Using the previous modules you can already use the functions CDFsphericity[p_,n_,x_,dist_:1] and PDFsphericity[p_,n_,x_,dist_:1] available in this file for the pdf and cdf of the near-exact distributions for the product of Generalized Gamma distributions *)
(* 5. Please see the instruction file for more details *)

CDFsphericity[p_,n_,x_,dist_:1]:=Module[{a,b,bstar,cstar,rj,m,r,vec1,vec2,vec3,vec4,vec5,vec6},a=Table[(n+1)/2-j/2,{j,2,p}];
Clear[vec1,vec2,vec3,vec4,vec5,vec6];
b=Table[(j-1)/p+(j-1)/2,{j,2,p}];
bstar=Floor[b];cstar=Rationalize[b-bstar,0];
rj=Flatten[{If[Mod[p,2]==0,Floor[p/2+1/2],Floor[p/4]],Table[Floor[(p-j)/2+1],{j,3,p}]}];
If[dist==1,mm1=SetPrecision[MomSLB[a+bstar,cstar,1],50];
mm2=SetPrecision[MomSLB[a+bstar,cstar,2],50];
m=mm1/(mm2-mm1^2);r=mm1^2/(mm2-mm1^2);
GNIGcdf[rj,r,a,m,x],If[dist==2,{vec1,vec2,vec3,vec4}={vec1,vec2,vec3,vec4}/.NSolve[Table[SetPrecision[MomSLB[a+bstar,cstar,h],500],{h,1,4}]==Table[Moment[MixtureDistribution[{vec4,1-vec4},{GammaDistribution[vec1,1/vec3],GammaDistribution[vec2,1/vec3]}],h],{h,1,4}],{vec1,vec2,vec3,vec4},Reals][[1]];
vec4*GNIGcdf[rj,vec1,a,vec3,x]+(1-vec4)*GNIGcdf[rj,vec2,a,vec3,x],{vec1,vec2,vec3,vec4,vec5,vec6}=Select[{vec1,vec2,vec3,vec4,vec5,vec6}/.NSolve[Table[SetPrecision[MomSLB[a+bstar,cstar,h],500],{h,1,6}]==Table[Moment[MixtureDistribution[{vec5,vec6,1-vec5-vec6},{GammaDistribution[vec1,1/vec4],GammaDistribution[vec2,1/vec4],GammaDistribution[vec3,1/vec4]}],h],{h,1,6}],{vec1,vec2,vec3,vec4,vec5,vec6},Reals],AllTrue[#>0&]][[1]];
vec5*GNIGcdf[rj,vec1,a,vec4,x]+vec6*GNIGcdf[rj,vec2,a,vec4,x]+(1-vec5-vec6)*GNIGcdf[rj,vec3,a,vec4,x]]]]


PDFsphericity[p_,n_,x_,dist_:1]:=Module[{a,b,bstar,cstar,rj,m,r,vec1,vec2,vec3,vec4,vec5,vec6},a=Table[(n+1)/2-j/2,{j,2,p}];
Clear[vec1,vec2,vec3,vec4,vec5,vec6];
b=Table[(j-1)/p+(j-1)/2,{j,2,p}];
bstar=Floor[b];cstar=Rationalize[b-bstar,0];
rj=Flatten[{If[Mod[p,2]==0,Floor[p/2+1/2],Floor[p/4]],Table[Floor[(p-j)/2+1],{j,3,p}]}];
If[dist==1,mm1=SetPrecision[MomSLB[a+bstar,cstar,1],50];
mm2=SetPrecision[MomSLB[a+bstar,cstar,2],50];
m=mm1/(mm2-mm1^2);r=mm1^2/(mm2-mm1^2);
GNIGpdf[rj,r,a,m,x],If[dist==2,{vec1,vec2,vec3,vec4}={vec1,vec2,vec3,vec4}/.NSolve[Table[SetPrecision[MomSLB[a+bstar,cstar,h],500],{h,1,4}]==Table[Moment[MixtureDistribution[{vec4,1-vec4},{GammaDistribution[vec1,1/vec3],GammaDistribution[vec2,1/vec3]}],h],{h,1,4}],{vec1,vec2,vec3,vec4},Reals][[1]];
vec4*GNIGpdf[rj,vec1,a,vec3,x]+(1-vec4)*GNIGpdf[rj,vec2,a,vec3,x],{vec1,vec2,vec3,vec4,vec5,vec6}=Select[{vec1,vec2,vec3,vec4,vec5,vec6}/.NSolve[Table[SetPrecision[MomSLB[a+bstar,cstar,h],500],{h,1,6}]==Table[Moment[MixtureDistribution[{vec5,vec6,1-vec5-vec6},{GammaDistribution[vec1,1/vec4],GammaDistribution[vec2,1/vec4],GammaDistribution[vec3,1/vec4]}],h],{h,1,6}],{vec1,vec2,vec3,vec4,vec5,vec6},Reals],AllTrue[#>0&]][[1]];
vec5*GNIGpdf[rj,vec1,a,vec4,x]+vec6*GNIGpdf[rj,vec2,a,vec4,x]+(1-vec5-vec6)*GNIGpdf[rj,vec3,a,vec4,x]]]]

MomSLB[a_,b_,h_]:=D[MGFSLB[a,b,t],{t,h}]/.t->0

MGFSLB[a_,b_,t_]:=Product[Gamma[a[[j]]-t]/Gamma[a[[j]]]*Gamma[a[[j]]+b[[j]]]/Gamma[a[[j]]+b[[j]]-t],{j,1,Length[a]}]

SetPrecision[1-CDFsphericity[5,50,Rationalize[-Log[0.6109257783234],0],1],20]

SetPrecision[1-CDFsphericity[5,50,Rationalize[-Log[0.6109257783234],0],2],20]

SetPrecision[1-CDFsphericity[5,50,Rationalize[-Log[0.6109257783234],0],3],20]