Minimizing NExpectation for a custom distribution in Mathematica

Solution 1:

As far as I see, the problem is (as you already wrote), that MeanResidualLife takes a long time to compute, even for a single evaluation. Now, the FindMinimum or similar functions try to find a minimum to the function. Finding a minimum requires either to set the first derivative of the function zero and solve for a solution. Since your function is quite complicated (and probably not differentiable), the second possibility is to do a numerical minimization, which requires many evaluations of your function. Ergo, it is very very slow.

I'd suggest to try it without Mathematica magic.

First let's see what the MeanResidualLife is, as you defined it. NExpectation or Expectation compute the expected value. For the expected value, we only need the PDF of your distribution. Let's extract it from your definition above into simple functions:

pdf[a_, b_, m_, s_, x_] := (1/(2*(a + b)))*a*b*
    (E^(a*(m + (a*s^2)/2 - x))*Erfc[(m + a*s^2 - x)/(Sqrt[2]*s)] + 
    E^(b*(-m + (b*s^2)/2 + x))*Erfc[(-m + b*s^2 + x)/(Sqrt[2]*s)])
pdf2[a_, b_, m_, s_, x_] := pdf[a, b, m, s, Log[x]]/x;

If we plot pdf2 it looks exactly as your Plot

Plot[pdf2[3.77, 1.34, -2.65, 0.40, x], {x, 0, .3}]

Now to the expected value. If I understand it correctly we have to integrate x * pdf[x] from -inf to +inf for a normal expected value.

x * pdf[x] looks like

Plot[pdf2[3.77, 1.34, -2.65, 0.40, x]*x, {x, 0, .3}, PlotRange -> All]

and the expected value is

NIntegrate[pdf2[3.77, 1.34, -2.65, 0.40, x]*x, {x, 0, \[Infinity]}]
Out= 0.0596504

But since you want the expected value between a start and +inf we need to integrate in this range, and since the PDF then no longer integrates to 1 in this smaller interval, I guess we have to normalize the result be dividing by the integral of the PDF in this range. So my guess for the left-bound expected value is

expVal[start_] := 
    NIntegrate[pdf2[3.77, 1.34, -2.65, 0.40, x]*x, {x, start, \[Infinity]}]/
    NIntegrate[pdf2[3.77, 1.34, -2.65, 0.40, x], {x, start, \[Infinity]}]

And for the MeanResidualLife you subtract start from it, giving

MRL[start_] := expVal[start] - start

Which plots as

Plot[MRL[start], {start, 0, 0.3}, PlotRange -> {0, All}]

Looks plausible, but I'm no expert. So finally we want to minimize it, i.e. find the start for which this function is a local minimum. The minimum seems to be around 0.05, but let's find a more exact value starting from that guess

FindMinimum[MRL[start], {start, 0.05}]

and after some errors (your function is not defined below 0, so I guess the minimizer pokes a little in that forbidden region) we get

{0.0418137, {start -> 0.0584312}}

So the optimum should be at start = 0.0584312 with a mean residual life of 0.0418137.

I don't know if this is correct, but it seems plausible.