0
$\begingroup$

I would like to estimate the parameters of the following distribution:

fSurvivalGompertzPDF[\[Alpha]_, \[Beta]_] := ProbabilityDistribution[(1/((
 E^(\[Alpha]/\[Beta]) Gamma[0, \[Alpha]/\[Beta]])/\[Beta])
  E^(((1 - E^(t \[Beta])) \[Alpha])/\[Beta])), {t, 0, \[Infinity]}]

If I create some data using this distribution:

data=RandomVariate[fSurvivalGompertzPDF[0.16, 0.65], {1000}]

Then try to estimate its parameters using:

FindDistributionParameters[RandomVariate[data, fSurvivalGompertzPDF[[Alpha]1, \[Beta]1]]

It hangs for an eternity. I have tried reducing the working precision, but this doesn't help. I have also tried to find its moments (to construct method of moments estimators for the parameters), but even $Mean$ and $Variance$ appear to be getting stuck.

I also tried doing the integrals manually, and I get the following:

  mean=MeijerG[{{}, {1, 1}}, {{0, 0, 0}, {}}, \[Alpha]/\[Beta]]/(\[Beta] Gamma[0, \[Alpha]/\[Beta]])
 variance=(2 MeijerG[{{}, {1, 1, 1}}, {{0, 0, 0, 
 0}, {}}, \[Alpha]/\[Beta]])/(\[Beta]^2 Gamma[0, \[Alpha]/\[Beta]]) - (MeijerG[{{}, {1, 1}}, {{0, 0, 0}, {}}, \[Alpha]/\[Beta]]/(\[Beta] Gamma[
   0, \[Alpha]/\[Beta]]))^2

Which are pretty foreboding, given that they have the Meijer G function! Predictably, trying to use the above moments to solve for $\alpha$ and $\beta$ with $NSolve$ doesn't work - it just returns the whole thing unevaluated

Not sure what else to try. I know constructed the log-likelihood manually will not produce answers in a reasonable running time (see a previous question by me.)

Does anyone have any ideas here?

Best,

Ben

$\endgroup$
1
$\begingroup$

Your fSurvivalGompertzPDF is a distribution not a PDF so I have renamed it to fSurvivalGompertzDist. The ProbabilityDistribution should include the distribution parameter assumptions to facilitate use of other built-in functions.

assume = {α > 0, β > 0};

fSurvivalGompertzDist[α_, β_] :=
 ProbabilityDistribution[
  (1/((E^(α/β) Gamma[0, α/β])/β) *
    E^(((1 - E^(t β)) α)/β)), {t, 0, ∞},
  Assumptions -> assume]

The PDF is then

pdf[t_] = PDF[fSurvivalGompertzDist[α, β], t]

(*  Piecewise[{{(E^(-(α/β) + ((1 - 
              E^(t*β))*α)/β)*
            β)/Gamma[0, α/β], t > 0}}, 0]  *)

Since the definition of the distribution includes the distribution parameter assumptions, Mean and Variance and CDF will use these assumptions

μ = Mean[fSurvivalGompertzDist[α, β]]

(*  MeijerG[{{}, {1, 1}}, {{0, 0, 0}, {}}, α/β]/(β Gamma[
  0, α/β])  *)

var = Variance[fSurvivalGompertzDist[α, β]]

(*  (1/(β^2 Gamma[
  0, α/β]^2))(-MeijerG[{{}, {1, 1}}, {{0, 0, 
      0}, {}}, α/β]^2 + 
  2 Gamma[0, α/β] MeijerG[{{}, {1, 1, 1}}, {{0, 0, 0, 
      0}, {}}, α/β])  *)

The CDF is

CDF[fSurvivalGompertzDist[α, β], t]

(*  Piecewise[{{(-ExpIntegralEi[-(α/β)] + 
            ExpIntegralEi[-((E^(t*β)*α)/β)])/
         Gamma[0, α/β], t > 0}}, 0]  *)

However, this can be simplified

Clear[cdf]

cdf[t_] = Assuming[assume, 
  CDF[fSurvivalGompertzDist[α, β], t] // FunctionExpand // 
   Simplify]

(*  Piecewise[{{1 - ExpIntegralEi[-((E^(t*β)*α)/β)]/
           ExpIntegralEi[-(α/β)], t > 0}}, 0]  *)

Verifying that the PDF is the derivative of the CDF

Assuming[And @@ assume && t > 0, pdf[t] == cdf'[t] // FullSimplify]

(*  True  *)

However, RandomVariate is not defined for this user-defined distribution

RandomVariate[fSurvivalGompertzDist[0.16, 0.65]]

(*  RandomVariate[
 ProbabilityDistribution[
  0.481049 E^(0.246154 (1 - E^(0.65 \[FormalX]))), {\[FormalX], 
   0, ∞}, Assumptions -> {α > 0, β > 0}]]  *)

Consequently, your data is not a list of positive real numbers and you must define your own random variate generator.

You can use InverseCDF to generate the random variates

SetAttributes[randVar, HoldFirst]

randVar[fSurvivalGompertzDist[α_, β_], {n_}] :=
 Table[InverseCDF[fSurvivalGompertzDist[α, β], 
   RandomReal[]], {n}]

randVar[fSurvivalGompertzDist[0.16, 0.65], {5}]

(*  {1.22391, 1.16153, 1.89787, 0.658346, 3.20399}  *)

Generating 1000 data values takes considerable time.

(data = randVar[fSurvivalGompertzDist[0.16, 0.65], {1000}];) // AbsoluteTiming

(*  {381.949, Null}  *)

With a valid set of data, FindDistributionParameters performs as expected.

FindDistributionParameters[data, fSurvivalGompertzDist[α, β]]

(*  {α -> 0.140665, β -> 0.670161}  *)
$\endgroup$
  • $\begingroup$ But RandomVariate works for fSurvivalGomperzDist when the Assumptions are not included in its definition. I'm sure this something else about Mathematica that I don't understand. Why would being more specific remove a feature? $\endgroup$ – JimB Feb 22 '16 at 3:37
  • 1
    $\begingroup$ @JimBaldwin - presumably, the definition of RandomVariate isn't flexible enough to cover the case when ProbabilityDistribution has the option Assumptions (or perhaps other options). Consequently, there is a trade-off when using the option. The specific intended use would determine which approach represents the path of least resistance. Or there may be another workaround. $\endgroup$ – Bob Hanlon Feb 22 '16 at 4:01
5
$\begingroup$

I don't know why maximizing the log likelihood didn't work for you. Below is an example with 1,000 random samples that takes just a few seconds:

(* Generate some data *)
data = RandomVariate[fSurvivalGompertzPDF[1, 1], 1000];

(* Construct log likelihood *)
logL[α_, β_] := 
 LogLikelihood[fSurvivalGompertzPDF[α, β], data]

(* Find values that maximize the log likelihood *)
FindMaximum[{logL[α, β], α > 0 && β > 0}, {{α, 1}, {β, 1}}]
(* {-141.24890424996966,{α->0.8324136054901223,β->1.2610928114287494}} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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