(*****************************************************************************************************************)
(* 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 nearexactPDF[r_,lambda_,beta_,gamma_,m_,w_] and nearexactCDF[r_,lambda_,beta_,gamma_,m_,w_] 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 *)
	

(* Mixture of Gamma *)

FCSMG[p_,rrr_,la_,ww_,t_]:=Module[{},
c=Length[p];
pp=Sum[p[[j]],{j,1,c}];
Sum[p[[j]]*la[[j]]^rrr[[j]]*(la[[j]]-I*t)^(-rrr[[j]])*Exp[ww[[j]]*I*t],{j,1,c}]+(1-pp)*la[[1+c]]^rrr[[1+c]]*(la[[1+c]]-I*t)^(-rrr[[1+c]])*Exp[ww[[1+c]]*I*t]]

MomSMG[p_,rrr_,la_,ww_,h_]:=I^(-h)*D[FCSMG[p,rrr,la,ww,t],{t,h}]/.t->0

(* PDF FOR THE NEAR-EXACT DISTRIBUTIONS OF THE PRODUCT OF INDEPENDENT GENERALIZED GAMMA RANDOM VARIABLES *)

nearexactPDF[r_,lambda_,beta_,gamma_,m_,w_]:=Module[{mom,mom1,mom2,v,vs,n,isc,rr,shift,Lambda},
ClearAll[omega,rho,upsilon,pe,ppp,pppp];
pe=Table[Unique["pe"],{m}];
mom=Table[SetPrecision[I^(-h)*D[Product[Gamma[r[[j]]+gamma-I*t/beta[[j]]]/Gamma[r[[j]]+gamma],{j,1,Length[r]}],{t,h}]/.t->0,500],{h,1,m}];
mom1=Table[SetPrecision[I^(-h)*D[omega^rho*(omega-I*t)^(-rho)*Exp[I*t*upsilon],{t,h}]/.t->0,500],{h,1,3}];
{rho,omega,upsilon}={rho,omega,upsilon}/.Flatten[NSolve[{mom[[1]]==mom1[[1]],mom[[2]]==mom1[[2]],mom[[3]]==mom1[[3]]},{rho,omega,upsilon}]];
mom2=SetPrecision[Table[MomSMG[pe,Table[rho+i,{i,0,m}],Table[omega,{i,1,m+1}],Table[upsilon,{i,1,m+1}],h],{h,1,m}],50];
ppp=Flatten[pe/.
NSolve[Table[mom[[h]]==mom2[[h]],{h,1,m}],pe]];
pppp=Flatten[{ppp,{1-Sum[ppp[[k]],{k,1,m}]}}];
shift=upsilon+Log[Product[lambda[[j]]^(1/beta[[j]]),{j,1,Length[r]}]];
v=Flatten[{Table[Table[((r[[j]]+k)*beta[[j]]),{k,0,gamma-1}],{j,1,Length[r]}]}];
vs=Sort[v];n=Length[v];Lambda={vs[[1]]};rr={1};isc=1;
Do[If[vs[[i]]==vs[[i-1]],{rr[[isc]]=rr[[isc]]+1},{isc=isc+1,Lambda=Append[Lambda,vs[[i]]],rr=Append[rr,1]}],{i,2,n}];
Sum[pppp[[k]]*(GNIGpdf[rr,rho+(k-1),Lambda,omega,-Log[w]-shift])*1/w,{k,1,m+1}]
]

(* CDF FOR THE NEAR-EXACT DISTRIBUTIONS OF THE PRODUCT OF INDEPENDENT GENERALIZED GAMMA RANDOM VARIABLES *)

nearexactCDF[r_,lambda_,beta_,gamma_,m_,w_]:=Module[{mom,mom1,mom2,v,vs,n,isc,rr,shift,Lambda},
ClearAll[omega,rho,upsilon,pe,ppp,pppp];
pe=Table[Unique["pe"],{m}];
mom=Table[SetPrecision[I^(-h)*D[Product[Gamma[r[[j]]+gamma-I*t/beta[[j]]]/Gamma[r[[j]]+gamma],{j,1,Length[r]}],{t,h}]/.t->0,500],{h,1,m}];
mom1=Table[SetPrecision[I^(-h)*D[omega^rho*(omega-I*t)^(-rho)*Exp[I*t*upsilon],{t,h}]/.t->0,500],{h,1,3}];
{rho,omega,upsilon}={rho,omega,upsilon}/.Flatten[NSolve[{mom[[1]]==mom1[[1]],mom[[2]]==mom1[[2]],mom[[3]]==mom1[[3]]},{rho,omega,upsilon}]];
mom2=SetPrecision[Table[MomSMG[pe,Table[rho+i,{i,0,m}],Table[omega,{i,1,m+1}],Table[upsilon,{i,1,m+1}],h],{h,1,m}],50];
ppp=Flatten[pe/.NSolve[Table[mom[[h]]==mom2[[h]],{h,1,m}],pe]];
pppp=Flatten[{ppp,{1-Sum[ppp[[k]],{k,1,m}]}}];
shift=upsilon+Log[Product[lambda[[j]]^(1/beta[[j]]),{j,1,Length[r]}]];
v=Flatten[{Table[Table[((r[[j]]+k)*beta[[j]]),{k,0,gamma-1}],{j,1,Length[r]}]}];
vs=Sort[v];n=Length[v];Lambda={vs[[1]]};rr={1};isc=1;
Do[If[vs[[i]]==vs[[i-1]],{rr[[isc]]=rr[[isc]]+1},{isc=isc+1,Lambda=Append[Lambda,vs[[i]]],rr=Append[rr,1]}],{i,2,n}];
1-Sum[pppp[[k]]*(GNIGcdf[rr,rho+(k-1),Lambda,omega,-Log[w]-shift]),{k,1,m+1}]]

(* Example *)

Clear[r,lambda,beta,gamma,m]
r={2,3,5};
lambda={1,2,10};
beta={5,6,7};
gamma=4;
m=4;

(* Quantiles *)
q1=w/.FindRoot[nearexactCDF[r,lambda,beta,gamma,m,w]==90/100,{w,1.1},WorkingPrecision->30]
q2=w/.FindRoot[nearexactCDF[r,lambda,beta,gamma,m,w]==95/100,{w,1.5},WorkingPrecision->30]
q3=w/.FindRoot[nearexactCDF[r,lambda,beta,gamma,m,w]==99/100,{w,1.8},WorkingPrecision->30]

(* Using the CDF *)

SetPrecision[nearexactCDF[r,lambda,beta,gamma,m,q1],15]
SetPrecision[nearexactCDF[r,lambda,beta,gamma,m,q2],15]
SetPrecision[nearexactCDF[r,lambda,beta,gamma,m,q3],15]

(* Using the PDF *)

SetPrecision[nearexactPDF[r,lambda,beta,gamma,m,q1],15]
SetPrecision[nearexactPDF[r,lambda,beta,gamma,m,q2],15]
SetPrecision[nearexactPDF[r,lambda,beta,gamma,m,q3],15]

fig=Plot[nearexactPDF[r,lambda,beta,gamma,m,w],{w,.00001,2}]

hist=Histogram[Product[RandomVariate[GammaDistribution[r[[j]],1/lambda[[j]]],800000]^(1/beta[[j]]),{j,1,Length[r]}],Automatic,"ProbabilityDensity"]

Show[fig,hist ]