Suppose the i.i.d. sample Y is fitted with models A, B, and C.
In order to determine which model’s prediction is likely to be closer to the true (but unknown) distribution that generated the original sample Y, one of the standard prescriptions is to use the AIC. Select the model with the smallest AIC.
Ref. Akaike information criterion - Wikipedia
In the following example, models A, B, and C are respectively Gamma, Weibull, and LogNormal. I have experimented with three true distributions Gamma(5, 2), Weibull(2.4, 11.3), and LogNormal(2.2, 0.47). These are similar each others. Please look at the graph in the following example.
Results of a Monte Carlo simulation (10^4 iterations) with samples Y generated by Gamma(5, 2) of size 100:
Gamma => 0.6692
Weibull => 0.1409
LogNormal => 0.1899
The percentage of the correct answer, Gamma, being selected is about 67 percent.
In the cases of Weibull(2.4, 11.3) and LogNormal(2.2, 0.47), the percentages of the correct answers are about 86 percent and about 81 percent, respectively.
For details, see the following code and results.
using Distributions, StatsPlots
function aic(model, Y)
mle = fit_mle(model, Y)
-2loglikelihood(mle, Y) + 2length(params(mle))
end
function simulate_model_selections(models, truedist, samplesize; niters = 10^4)
selectedmodel = Vector{Int}(undef, niters)
Threads.@threads for i in 1:niters
Y = rand(truedist, samplesize)
selectedmodel[i] = argmin(aic.(models, Ref(Y)))
end
nselected = zeros(Int, 3)
for i in 1:niters
nselected[selectedmodel[i]] += 1
end
[model => nselected[i]/niters for (i, model) in enumerate(models)]
end
models = (Gamma, Weibull, LogNormal)
a, b = 0, 40
plot(Gamma(5, 2), a, b; label="Gamma(5, 2)")
plot!(Weibull(2.4, 11.3), a, b; label="Weibull(2.4, 11.3)", ls=:dash)
plot!(LogNormal(2.2, 0.47), a, b; label="LogNormal(2.2, 0.47)", ls=:dashdot)

simulate_model_selections(models, Gamma(5, 2), 100; niters = 10^4)
3-element Vector{Pair{UnionAll, Float64}}:
Gamma => 0.6692
Weibull => 0.1409
LogNormal => 0.1899
simulate_model_selections(models, Weibull(2.4, 11.3), 100; niters = 10^4)
3-element Vector{Pair{UnionAll, Float64}}:
Gamma => 0.1427
Weibull => 0.856
LogNormal => 0.0013
simulate_model_selections(models, LogNormal(2.2, 0.47), 100; niters = 10^4)
3-element Vector{Pair{UnionAll, Float64}}:
Gamma => 0.1914
Weibull => 0.0034
LogNormal => 0.8052
Pluto notebook: Simulation of AIC model selection - Pluto notebook.jl
You can execute in Pluto by copy-and-pasting the URL into “Open from file:” input form and click “Open” button.

Edit 1: I’m sorry. The original code works with Julia ≥ v1.7.0-beta. The edited code above works also on Julia v1.6.1.
Edit 2: Link to Pluto notebook.