2
$\begingroup$

I write the density function of cw distribution and add data to calculate parameters maximum likelihood but the results of Mathematica no like the researchcode Mathematica for calculate parameters

dist = ProbabilityDistribution[(1 - Exp[λ*(1 - Exp[1])])^(-1) *λ*β*(γ/α)*(x/α)^(γ - 1)*
      (1 - Exp[-(x/α)^γ])^(β - 1)*Exp[λ*(1 - Exp[-(x/α)^γ])^β - (x/α)^γ + (1 - Exp[-(x/α)^γ])], {x, 1, ∞}, 
      Assumptions -> {λ > 0, β > 0, α > 0, γ > 0}]
data = { 70,  90,  96,  97,  99, 103, 104, 104, 105, 107, 108, 108, 108, 109,
        109 ,112, 112 ,113 ,114 ,114 ,114, 116 ,119 ,120, 120, 120, 121, 121,
        123 ,124, 124 ,124, 124, 124 ,128, 128 ,129 ,129 ,130, 130, 130, 131,
        131, 131 ,131 ,131 ,132, 132, 132 ,133, 134, 134, 134 ,134 ,136 ,136,
        137, 138 ,138, 138, 139, 139, 141, 141 ,142 ,142 ,142, 142, 142 ,142,
        142, 142, 144 ,144 ,145 ,146 ,148 ,148, 149, 151 ,151 ,152 ,155, 156,
        155 ,156 ,157 ,157 ,157 ,157, 158 ,159, 162, 163 ,163, 164 ,166 ,166,
        168, 170 ,174, 201 ,212}

I want the values of parameters as below Could you help me to find codethe values of parameters

$\endgroup$
0

3 Answers 3

3
$\begingroup$

You have an error in your pdf. It should be the following:

pdf = (E^((1 - E^-(x/α)^γ)^β - (x/α)^γ + (1 - E^(1 - E^-(x/α)^γ)^β) λ) (1 - E^-(x/α)^γ)^(-1 + β)*
  (x/α)^(-1 + γ) β γ λ)/((1 - E^(λ - E λ)) α)

(See reference.) Also, the limits of integration are from 0 to $\infty$ rather than from 1 to $\infty$.

To obtain the maximum likelihood estimates of the parameters:

dist = ProbabilityDistribution[pdf, {x, 0, ∞}, 
   Assumptions -> {λ > 0, α > 0, β > 0, η > 0}];
data = {70, 90, 96, 97, 99, 103, 104, 104, 105, 107, 108, 108, 108, 
   109, 109, 112, 112, 113, 114, 114, 114, 116, 119, 120, 120, 120, 
   121, 121, 123, 124, 124, 124, 124, 124, 128, 128, 129, 129, 130, 
   130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 133, 134, 134, 
   134, 134, 136, 136, 137, 138, 138, 138, 139, 139, 141, 141, 142, 
   142, 142, 142, 142, 142, 142, 142, 144, 144, 145, 146, 148, 148, 
   149, 151, 151, 152, 155, 156, 155, 156, 157, 157, 157, 157, 158, 
   159, 162, 163, 163, 164, 166, 166, 168, 170, 174, 201, 212};
mle = FindDistributionParameters[data, dist, {{α, 133}, {β, 3.45}, {γ, 3.3}, {λ, 3}},
  ParameterEstimator -> "MaximumLikelihood"]
(* {α -> 133.646, β -> 3.45368, γ -> 3.37814, λ -> 2.91285} *)

Show[SmoothHistogram[data, 
  PlotLegends -> LineLegend[{Blue, Red}, {"Smoothed histogram", "PDF with MLE"}]],
  Plot[pdf /. mle, {x, 0, 250}, PlotRange -> All, PlotStyle -> Red]]

Smoothed histogram and pdf with mle

I gave up waiting for the "MethodOfMoments" estimates.

You might ask why are these estimates so different from what you expected? In part that's because some of the estimators are highly correlated with each other:

(* Log of the likelihood *)
logL = LogLikelihood[dist, data];

(* Parameter covariance matrix *)
(covMat = -Inverse[(D[logL, {{α, β, γ, λ}, 2}]) /. mle]) // MatrixForm

Parameter covariance matrix

(* Parameter correlation matrix *)
(corMat = Table[covMat[[i, j]]/(Sqrt[covMat[[i, i]] covMat[[j, j]]]), 
      {i, 4}, {j, 4}]) // MatrixForm

Parameter correlation matrix

We see the the estimators for $\alpha$, $\beta$, and $\gamma$ are highly correlated with each other and the standard errors for the parameters are large compared to the estimates:

Sqrt[Diagonal[covMat]]
(* {23.6787, 2.6356, 1.48478, 1.43782} *)

And maybe the data just doesn't warrant something more complex than a simple Weibull distribution.

$\endgroup$
7
  • $\begingroup$ Why you selected {α, 133}, {β, 3.45}, {γ, 3.3}, {λ, 3} in finddistributionparameter $\endgroup$
    – A Day
    Mar 1, 2021 at 19:16
  • $\begingroup$ One needs to choose some starting values because without specifying starting values FindDistributionParameters returns a negative value for the estimate of $\gamma$. Choosing {{\[Alpha], 1}, {\[Beta], 1}, {\[Gamma], 1}, {\[Lambda], 1}} also works. $\endgroup$
    – JimB
    Mar 1, 2021 at 19:58
  • $\begingroup$ Is this affect on log likelihood and aic $\endgroup$
    – A Day
    Mar 1, 2021 at 20:40
  • $\begingroup$ How I compute best estimators $\endgroup$
    – A Day
    Mar 1, 2021 at 20:40
  • $\begingroup$ If the starting values end up at the global maximum of the log likelihood function by definition there is no effect on the log likelihood and aic. I don't know if there are other local maxima that could cause trouble. There's probably more "trouble" introduced by your data being rounded to the nearest integer. Is there a reason for rounding? $\endgroup$
    – JimB
    Mar 1, 2021 at 21:18
2
$\begingroup$

An extended comment:

  1. The expression in your ProbabilityDistribution is very much not a PDF, i.e. it does not integrate to $1$ even with your best-fit parameters:

    NIntegrate[ PDF[dist /. {λ -> 6.775, β -> 35.209, α -> 49.29, γ -> 1.019}][x], {x, 0, 550}]
    (* Out: 2360.9 *)
    

    Adding Method -> "Normalize" to the ProbabilityDistribution expression makes it extremely slow, to the point that one can't obtain an answer in a reasonable time.

  2. Just by inspection, though, the parameter values in your table do NOT seem to be a good fit to your data. First, here is the histogram for your data:

    Histogram[data, Automatic, "PDF", PlotRange -> {{0, 550}, Automatic}]
    

    histogram from experimental data

    ... and here is the plot of the PDF of your ProbabilityDistribution when your best-fit values are substituted in, plotted on the same horizontal scale:

    Plot[PDF[dist /. {λ -> 6.775, β -> 35.209, α -> 49.29, γ -> 1.019}][x], {x, 1, 550}]
    

    plot of PDF with "best-fit" values plugged in

    That plot looks like a pretty poor match to the histogram, on both the horizontal scale (the position of the maximum is off) and on the vertical scale (the PDF is not appropriately normalized).

$\endgroup$
1
$\begingroup$

Given that you've rounded the data to the nearest integer, pretending that the rounding doesn't matter might or might not result in not-so-hot results. In such cases the rounding can be explicitly acknowledged using the cdf rather than the pdf.

In short the likelihood after rounding to the nearest integer is the following:

$$L=\prod_{i=1}^n (F(x_i+1/2)-F(x_i-1/2))$$

where $F$ is the cdf (cumulative distribution function) and the log of the likelihood is

$$\log{(L)}=\sum_{i=1}\log{(F(x_i+1/2)-F(x_i-1/2))}$$

The maximum likelihood estimates are the values that maximize $\log{(L)}$:

(* data *)
data = {70, 90, 96, 97, 99, 103, 104, 104, 105, 107, 108, 108, 108, 
   109, 109, 112, 112, 113, 114, 114, 114, 116, 119, 120, 120, 120, 
   121, 121, 123, 124, 124, 124, 124, 124, 128, 128, 129, 129, 130, 
   130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 133, 134, 134, 
   134, 134, 136, 136, 137, 138, 138, 138, 139, 139, 141, 141, 142, 
   142, 142, 142, 142, 142, 142, 142, 144, 144, 145, 146, 148, 148, 
   149, 151, 151, 152, 155, 156, 155, 156, 157, 157, 157, 157, 158, 
   159, 162, 163, 163, 164, 166, 166, 168, 170, 174, 201, 212};

(* cdf *)
cdf[x_, α_, β_, γ_, λ_] := (1 - Exp[λ (1 - Exp[(1 - Exp[-(x/α)^γ])^β])])/(1 - E^(λ - E λ))

(* Log of the likelihood *)
logL = Total[Log[cdf[# + 0.5, α, β, γ, λ] - cdf[# - 0.5, α, β, γ, λ]] & /@ data];

(* Maximum likelihood estimates *)
sol = FindMaximum[{logL, α > 0 && β > 0 && γ > 0 && λ > 0}, 
  {{α, 133.65}, {β, 3.45}, {γ, 3.38}, {λ, 2.918}}]
(* {-463.416, {α -> 133.636, β -> 3.45469, γ -> 3.37792, λ -> 2.91281}} *)

(* Covariance matrix *)
(covMat = -Inverse[(D[logL, {{α, β, γ, λ}, 2}]) /. sol[[2]]]) // MatrixForm

Covariance matrix

(* Parameter standard errors *)
Sqrt[Diagonal[covMat]]
(* {23.6844, 2.63705, 1.48493, 1.43788} *)

(* Correlation matrix *)
(corMat = Table[covMat[[i, j]]/(Sqrt[covMat[[i, i]] covMat[[j, j]]]), 
  {i, 4}, {j, 4}]) // MatrixForm

Correlation matrix

(* AIC *)
aic = -2 sol[[1]] + 2*4
(* 934.8312477645201 *)

In this case accounting for the rounding doesn't change much.

$\endgroup$
2
  • $\begingroup$ But the parameters above I wrote gives aic=859.704 and ks=.092 and I want the code give me this results to compare with my distribution $\endgroup$
    – A Day
    Mar 2, 2021 at 6:17
  • $\begingroup$ It is pdf online called (Chen-G class of distributions) $\endgroup$
    – A Day
    Mar 2, 2021 at 6:24

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.