Remove["@"]

(*Generate data set of size n from Gamma distribution with parameters α0 and λ0*)
n=100;
α0=2.1;
λ0=1.41;
data=RandomVariate[GammaDistribution[α0,1/λ0],n];
Histogram[data]
μ1=Mean[data];
μ2=Mean[data^2];
(Subscript[(OverHat[α]), MM])=μ1^2/(μ2-μ1^2);
(Subscript[(OverHat[λ]), MM])=μ1/(μ2-μ1^2);
(Subscript[(OverHat[α]), MLE])=α/.FindRoot[Log[μ1]-Mean[Log[data]]==Log[α]-D[Gamma[α],α]/Gamma[α],{α,(Subscript[(OverHat[α]), MM])}];
(Subscript[(OverHat[λ]), MLE])=(Subscript[(OverHat[α]), MLE])/μ1;
Print["α0=",α0," λ0=",λ0];
Print["Methods of moments: (Subscript[(OverHat[α]), MM])=",(Subscript[(OverHat[α]), MM])," (Subscript[(OverHat[λ]), MM])=",(Subscript[(OverHat[λ]), MM])];
Print["Maximum likelihood: (Subscript[(OverHat[α]), MLE])=",(Subscript[(OverHat[α]), MLE])," (Subscript[(OverHat[λ]), MLE])=",(Subscript[(OverHat[λ]), MLE])];

(*True distribution α0 and λ0 *)
(*Generate m times a dataset of size n based on Gamma[α0,1/λ0], estimate parameters in each simulation*)
m=100000;
For[i=1, i≤m,i++,
datasim[i]=RandomVariate[GammaDistribution[α0,1/λ0],n];
mu1[i]=Mean[datasim[i]];
mu2[i]=Mean[datasim[i]^2];
αsim[i]=mu1[i]^2/(mu2[i]-mu1[i]^2);
λsim[i]=mu1[i]/(mu2[i]-mu1[i]^2);
];
plotdistalphamm=Histogram[Table[αsim[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]
plotdistlambdamm=Histogram[Table[λsim[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]

(*Methods of moments: Generate m times a dataset of size n based on Gamma[(Subscript[(OverHat[α]), MM]),1/(Subscript[(OverHat[λ]), MM])], estimate parameters in each simulation*)
m=100000;
For[i=1, i≤m,i++,
datasim[i]=RandomVariate[GammaDistribution[(Subscript[(OverHat[α]), MM]),1/(Subscript[(OverHat[λ]), MM])],n];
mu1[i]=Mean[datasim[i]];
mu2[i]=Mean[datasim[i]^2];
αsim[i]=mu1[i]^2/(mu2[i]-mu1[i]^2);
λsim[i]=mu1[i]/(mu2[i]-mu1[i]^2);
];
plotbootstrapdistalphamm=Histogram[Table[αsim[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]
plotbootstrapdistlambdamm=Histogram[Table[λsim[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]

(*Maximum Likelihood: Generate m times a dataset of size n based on Gamma[(OverHat[α]),1/(Subscript[(OverHat[λ]), MLE])], estimate parameters in each simulation*)
m=100000;
For[i=1, i≤m,i++,
datasim[i]=RandomVariate[GammaDistribution[(Subscript[(OverHat[α]), MLE]),1/(Subscript[(OverHat[λ]), MLE])],n];
mu1[i]=Mean[datasim[i]];
(Subscript[αsim, MLE])[i]=α/.FindRoot[Log[mu1[i]]-Mean[Log[datasim[i]]]==Log[α]-D[Gamma[α],α]/Gamma[α],{α,(Subscript[(OverHat[α]), MLE])}];
(Subscript[λsim, MLE])[i]=(Subscript[αsim, MLE])[i]/mu1[i];
];
plotbootsrapdistalphamle=Histogram[Table[(Subscript[αsim, MLE])[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]
plotbootstrapdistlambdamle=Histogram[Table[(Subscript[λsim, MLE])[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]

(*Maximum Likelihood: Generate m times a dataset of size n based on Gamma[α0,1/λ0], estimate parameters in each simulation*)
m=100000;
For[i=1, i≤m,i++,
datasim[i]=RandomVariate[GammaDistribution[α0,1/λ0],n];
mu1[i]=Mean[datasim[i]];
(Subscript[αsim, MLE])[i]=α/.FindRoot[Log[mu1[i]]-Mean[Log[datasim[i]]]==Log[α]-D[Gamma[α],α]/Gamma[α],{α,(Subscript[(OverHat[α]), MLE])}];
(Subscript[λsim, MLE])[i]=(Subscript[αsim, MLE])[i]/mu1[i];
];
plotdistalphamle=Histogram[Table[(Subscript[αsim, MLE])[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]
plotdistlambdamle=Histogram[Table[(Subscript[λsim, MLE])[i],{i,1,m}],PlotRange->{{0,4},Automatic},AxesOrigin->{0,0}]

(*Comparison of the distribution of the estimate for α and λ with MLE and the methods of moments if the sample size is n*)
Grid[{{"","MLE","MM"},{"λ",plotdistlambdamle,plotdistlambdamm},{"α",plotdistalphamle,plotdistalphamm}}]

(*Conclusion: Sampling distribution MLE less dispersed than Methods of moments estimates. MLE is probably more precise than methods of moments*)