# How to plot the the maximum likelihood estimates?

I would like to find the code for the graphical representation of the maximum likelihood estimates. For that first, I find the estimate values using the the code given below

data = {6, 18, 29, 28, 47, 55, 40, 150, 129, 184, 236, 237, 336, 219, 612, 434, 648, 706, 838, 1129, 1421, 118, 116, 1393, 1540, 1941, 2175, 2278, 2824, 2803, 26};
f[ξ_, γ_, δ_] := ProbabilityDistribution[
(ξ*γ*δ*x^(-γ - 1)*E^(-δ*x^-γ)*E^(ξ*E^(-δ*x^-γ)))/(E^ξ - 1)
, {x, 0, ∞}
, Assumptions -> ξ > 0 && γ > 0 && δ > 0
]
PDF[f[ξ, γ, δ], x]
FindDistributionParameters[data, f[ξ, γ, δ]]


The output obtained is {γ -> 0.713216, δ -> 11.9749, ξ -> 2.9519}.

I want to plot the log-likelihood function of above estimates as same as the sample figure shown below.

• Can you present data or its part? Commented Dec 5, 2022 at 9:47
• The following is the data. data = {6, 18, 29, 28, 47, 55, 40, 150, 129, 184, 236, 237, 336, 219, 612, 434, 648, 706, 838, 1129, 1421, 118, 116, 1393, 1540, 1941, 2175, 2278, 2824, 2803, 26} Commented Dec 5, 2022 at 9:49
• You should mention that each of the curves are evaluated with the remaining parameters at their maximum likelihood estimate.
– JimB
Commented Dec 5, 2022 at 16:30

This mainly duplicates the other two answers but with an automatic way to obtain the plot ranges and the construction of log likelihood contour plots. (I have not made an efficient single function to produce the two plots. If there were more than 3 parameters, then writing a function would be advisable.)

The initial calculations:

data = {6, 18, 29, 28, 47, 55, 40, 150, 129, 184, 236, 237, 336, 219,
612, 434, 648, 706, 838, 1129, 1421, 118, 116, 1393, 1540, 1941,
2175, 2278, 2824, 2803, 26};
f[ξ_, γ_, δ_] := ProbabilityDistribution[(ξ*γ*δ*x^(-γ - 1)*
E^(-δ*x^-γ)*E^(ξ*E^(-δ*x^-γ)))/(E^ξ - 1), {x, 0, ∞},
Assumptions -> ξ > 0 && γ > 0 && δ > 0];
mle = FindDistributionParameters[data, f[ξ, γ, δ]]


Separate maximum likelihood estimates, get the log of the likelihood, and construct the parameter covariance matrix. The covariance matrix will be used to determine the plotting range for each parameter.

mleξ = ξ /. mle
mleγ = γ /. mle
mleδ = δ /. mle
logL = LogLikelihood[f[ξ, γ, δ], data];
covmat = -Inverse[(D[logL, {{ξ, γ, δ}, 2}]) /. mle];


Generate the individual likelihood plots

Show[Plot[logL /. {γ -> mleγ, δ -> mleδ},
{ξ, mleξ - 2 Sqrt[covmat[[1, 1]]], mleξ + 2 Sqrt[covmat[[1, 1]]]},
Frame -> True, FrameLabel -> {"ξ", "Log of likelihood"}] ,
ListPlot[{{mleξ, logL /. mle}}, PlotStyle -> {{PointSize[0.02], Red}}]]

Show[Plot[logL /. {ξ -> mleξ, δ -> mleδ},
{γ, mleγ - 2 Sqrt[covmat[[2, 2]]], mleγ + 2 Sqrt[covmat[[2, 2]]]},
Frame -> True, FrameLabel -> {"γ", "Log of likelihood"}],
ListPlot[{{mleγ, logL /. mle}}, PlotStyle -> {{PointSize[0.02], Red}}]]

Show[Plot[logL /. {ξ -> mleξ, γ -> mleγ},
{δ, mleδ - 2 Sqrt[covmat[[3, 3]]], mleδ + 2 Sqrt[covmat[[3, 3]]]},
Frame -> True, FrameLabel -> {"δ", "Log of likelihood"}],
ListPlot[{{mleδ, logL /. mle}}, PlotStyle -> {{PointSize[0.02], Red}}]]


Display the associated contour plots.

Show[ContourPlot[logL /. δ -> mleδ,
{ξ, mleξ - 2 Sqrt[covmat[[1, 1]]], mleξ + 2 Sqrt[covmat[[1, 1]]]},
{γ, mleγ - 2 Sqrt[covmat[[2, 2]]], mleγ + 2 Sqrt[covmat[[2, 2]]]},
Frame -> True, FrameLabel -> {"ξ", "γ"}, PlotLabel -> "Log of the likelihood",
Contours -> 20, ContourLabels -> True],
ListPlot[{{mleξ, mleγ}}, PlotStyle -> {{PointSize[0.02], Red}}]]

Show[ContourPlot[logL /. γ -> mleγ,
{ξ, mleξ - 2 Sqrt[covmat[[1, 1]]], mleξ + 2 Sqrt[covmat[[1, 1]]]},
{δ, mleγ - 2 Sqrt[covmat[[3, 3]]], mleδ + 2 Sqrt[covmat[[3, 3]]]},
Frame -> True, FrameLabel -> {"ξ", "δ"}, PlotLabel -> "Log of the likelihood",
Contours -> 20, ContourLabels -> True],
ListPlot[{{mleξ, mleδ}},
PlotStyle -> {{PointSize[0.02], Red}}]]

Show[ContourPlot[logL /. ξ -> mleξ,
{γ, mleγ - 2 Sqrt[covmat[[2, 2]]], mleγ + 2 Sqrt[covmat[[2, 2]]]},
{δ, mleδ - 2 Sqrt[covmat[[3, 3]]], mleδ + 2 Sqrt[covmat[[3, 3]]]},
Frame -> True, FrameLabel -> {"γ", "δ"}, PlotLabel -> "Log of the likelihood",
Contours -> 20, ContourLabels -> True],
ListPlot[{{mleγ, mleδ}}, PlotStyle -> {{PointSize[0.02], Red}}]]


First, let us reproduce your results for completeness.

data = {6, 18, 29, 28, 47, 55, 40, 150, 129, 184, 236, 237, 336, 219,
612, 434, 648, 706, 838, 1129, 1421, 118, 116, 1393, 1540, 1941,
2175, 2278, 2824, 2803, 26};
f[\[Xi]_, \[Gamma]_, \[Delta]_] := ProbabilityDistribution[(\[Xi]*\[Gamma]*\[Delta]*
x^(-\[Gamma] - 1)*E^(-\[Delta]*x^-\[Gamma])*E^(\[Xi]*E^(-\[Delta]*x^-\[Gamma])))/
(E^\[Xi] - 1),{x, 0, \[Infinity]},Assumptions-> \[Xi] > 0 && \[Gamma] > 0&&\[Delta] > 0];
PDF[f[\[Xi], \[Gamma], \[Delta]], x]


$$\begin{array}{cc} \{ & \begin{array}{cc} \frac{\gamma \delta \xi x^{-\gamma -1} e^{\xi e^{\delta \left(-x^{-\gamma }\right)}-\delta x^{-\gamma }}}{e^{\xi }-1} & x>0 \\ 0 & \text{True} \\ \end{array} \\ \end{array}$$

FindDistributionParameters[data, f[\[Xi], \[Gamma], \[Delta]]]


{\[Gamma] -> 0.713216, \[Delta] -> 11.9749, \[Xi] -> 2.9519}

To be sure, let us verify PDF by

NIntegrate[PDF[f[\[Xi], \[Gamma], \[Delta]], x] /.
{\[Xi] -> 1, \[Gamma] -> 1, \[Delta] -> 1}, {x, 0, Infinity}]


1.0

OK. Unfortunately,

LogLikelihood[(E^(-x^-\[Gamma] \[Delta] + E^(-x^-\[Gamma] \[Delta]) \[Xi])
x^(-1 - \[Gamma]) \[Gamma] \[Delta] \[Xi])/(-1 + E^\[Xi]), data]


returns the input. This is not any bug. Every command has its limitation and no result is better than an incorrect one. In view of it we have to do some work, writing down the loglikelyhood function up to its definition by

Dimensions[data]


{31}

Total[Table[Log[(E^(-x^-\[Gamma] \[Delta] + E^(-x^-\[Gamma] \[Delta]) \[Xi])
x^(-1 - \[Gamma]) \[Gamma] \[Delta] \[Xi])/(-1 + E^\[Xi])] /. x -> data[[j]], {j, 1, 31}]];

Plot[Total[ Table[Log[(E^(-x^-\[Gamma] \[Delta] +
E^(-x^-\[Gamma] \[Delta]) \[Xi]) x^(-1 - \[Gamma]) \[Gamma] \[Delta] \[Xi])/
(-1 + E^\[Xi])] /. x -> data[[j]], {j, 1,  31}]] /. {\[Delta] -> 11.9749, \[Xi] -> 2.9519}, {\[Gamma], 0.1,  1}]


I leave other plots on your own.

• Beautiful and informative.
– user88712
Commented Dec 5, 2022 at 10:31
• LogLikelihood[f[\[Xi], \[Gamma], \[Delta]], data] works (as it is the distribution rather than the pdf that is required for the first argument).
– JimB
Commented Dec 5, 2022 at 16:28
• @JimB: Thank you for your valuable comment. Commented Dec 5, 2022 at 17:55
Clear["Global*"]

data = {6, 18, 29, 28, 47, 55, 40, 150, 129, 184, 236, 237, 336, 219, 612,
434, 648, 706, 838, 1129, 1421, 118, 116, 1393, 1540, 1941, 2175, 2278,
2824, 2803, 26};

f[ξ_, γ_, δ_] :=
ProbabilityDistribution[(ξ*γ*δ*x^(-γ - 1)*
E^(-δ*x^-γ)*
E^(ξ*E^(-δ*x^-γ)))/(E^ξ - 1), {x, 0, ∞},
Assumptions -> ξ > 0 && γ > 0 && δ > 0]

param = FindDistributionParameters[data, f[ξ, γ, δ]]

(* {γ -> 0.713216, δ -> 11.9749, ξ -> 2.9519} *)


The first argument of the LogLikelihood must be the distribution not the PDF. The LogLikelihood is (only skeleton shown)

(llf = LogLikelihood[f[ξ, γ, δ], data]) // Short[#, 6] &


The Subsets of the parameters for plotting (fix two and plot one) are

ss = Reverse@Subsets[param, {2}]

(* {{δ -> 11.9749, ξ -> 2.9519}, {γ -> 0.713216, ξ ->
2.9519}, {γ -> 0.713216, δ -> 11.9749}} *)


The shown plot ranges of the parameters are

rng = {{γ, 0, 5}, {δ, 11, 13}, {ξ, 1.5, 5}};


Plotting,

Column[Plot[llf /. #[[1]], #[[2]],
Frame -> True,
ImageSize -> Medium,
FrameLabel -> {#[[2, 1]], LogLikelihood},
Epilog -> {Red, AbsolutePointSize[4],
Point[{#[[2, 1]], llf} /. param]}] & /@
Transpose[{ss, rng}]]
`