5
$\begingroup$

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.

enter image description here

$\endgroup$
3
  • $\begingroup$ Can you present data or its part? $\endgroup$
    – user64494
    Dec 5, 2022 at 9:47
  • $\begingroup$ 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} $\endgroup$
    – Deepthy Gs
    Dec 5, 2022 at 9:49
  • $\begingroup$ You should mention that each of the curves are evaluated with the remaining parameters at their maximum likelihood estimate. $\endgroup$
    – JimB
    Dec 5, 2022 at 16:30

3 Answers 3

7
$\begingroup$

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}}]]

Individual log likelihood plots

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}}]]

Log likelihood contour plots

$\endgroup$
6
$\begingroup$

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}]

enter image description here

I leave other plots on your own.

$\endgroup$
3
  • $\begingroup$ Beautiful and informative. $\endgroup$
    – user88712
    Dec 5, 2022 at 10:31
  • 5
    $\begingroup$ LogLikelihood[f[\[Xi], \[Gamma], \[Delta]], data] works (as it is the distribution rather than the pdf that is required for the first argument). $\endgroup$
    – JimB
    Dec 5, 2022 at 16:28
  • $\begingroup$ @JimB: Thank you for your valuable comment. $\endgroup$
    – user64494
    Dec 5, 2022 at 17:55
4
$\begingroup$
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] &

enter image description here

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}]]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.