Estimate parameters from distribution created with ProbabilityDistribution

Question

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

fSurvivalGompertzPDF[α_, β_] := ProbabilityDistribution[(1/((
 E^(α/β) Gamma[0, α/β])/β)
  E^(((1 - E^(t β)) α)/β)), {t, 0, ∞}]

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, β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}, {}}, α/β]/(β Gamma[0, α/β])
 variance=(2 MeijerG[{{}, {1, 1, 1}}, {{0, 0, 0, 
 0}, {}}, α/β])/(β^2 Gamma[0, α/β]) - (MeijerG[{{}, {1, 1}}, {{0, 0, 0}, {}}, α/β]/(β Gamma[
   0, α/β]))^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

Answer 1

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}  *)

Answer 2

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}} *)