# Estimate parameters from distribution created with ProbabilityDistribution

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

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

• 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? – JimB Feb 22 '16 at 3:37
• @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. – Bob Hanlon Feb 22 '16 at 4:01

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