title: Minimizing NExpectation for a custom distribution in Mathematica
body: <p>This relates to an earlier question from back in June:</p>
  
  <p><a href="https://stackoverflow.com/questions/6296455/calculating-expectation-for-a-custom-distribution-in-mathematica">Calculating expectation for a custom distribution in Mathematica</a></p>
  
  <p>I have a custom mixed distribution defined using a second custom distribution following along the lines discussed by <code>@Sasha</code> in a number of answers over the past year.</p>
  
  <p>Code defining the distributions follows:</p>
  
  <pre><code>nDist /: CharacteristicFunction[nDist[a_, b_, m_, s_], 
     t_] := (a b E^(I m t - (s^2 t^2)/2))/((I a + t) (-I b + t));
  nDist /: PDF[nDist[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)]); 
  nDist /: CDF[nDist[a_, b_, m_, s_], 
     x_] := ((1/(2*(a + b)))*((a + b)*E^(a*x)* 
          Erfc[(m - x)/(Sqrt[2]*s)] - 
         b*E^(a*m + (a^2*s^2)/2)*Erfc[(m + a*s^2 - x)/(Sqrt[2]*s)] + 
         a*E^((-b)*m + (b^2*s^2)/2 + a*x + b*x)*
          Erfc[(-m + b*s^2 + x)/(Sqrt[2]*s)]))/ E^(a*x);         
  
  nDist /: Quantile[nDist[a_, b_, m_, s_], p_] :=  
   Module[{x}, 
     x /. FindRoot[CDF[nDist[a, b, m, s], x] == #, {x, m}] &amp; /@ p] /; 
    VectorQ[p, 0 &lt; # &lt; 1 &amp;]
  nDist /: Quantile[nDist[a_, b_, m_, s_], p_] := 
   Module[{x}, x /. FindRoot[CDF[nDist[a, b, m, s], x] == p, {x, m}]] /;
     0 &lt; p &lt; 1
  nDist /: Quantile[nDist[a_, b_, m_, s_], p_] := -Infinity /; p == 0
  nDist /: Quantile[nDist[a_, b_, m_, s_], p_] := Infinity /; p == 1
  nDist /: Mean[nDist[a_, b_, m_, s_]] := 1/a - 1/b + m;
  nDist /: Variance[nDist[a_, b_, m_, s_]] := 1/a^2 + 1/b^2 + s^2;
  nDist /: StandardDeviation[ nDist[a_, b_, m_, s_]] := 
    Sqrt[ 1/a^2 + 1/b^2 + s^2];
  nDist /: DistributionDomain[nDist[a_, b_, m_, s_]] := 
   Interval[{0, Infinity}]
  nDist /: DistributionParameterQ[nDist[a_, b_, m_, s_]] := ! 
    TrueQ[Not[Element[{a, b, s, m}, Reals] &amp;&amp; a &gt; 0 &amp;&amp; b &gt; 0 &amp;&amp; s &gt; 0]]
  nDist /: DistributionParameterAssumptions[nDist[a_, b_, m_, s_]] := 
   Element[{a, b, s, m}, Reals] &amp;&amp; a &gt; 0 &amp;&amp; b &gt; 0 &amp;&amp; s &gt; 0
  nDist /: Random`DistributionVector[nDist[a_, b_, m_, s_], n_, prec_] :=
  
      RandomVariate[ExponentialDistribution[a], n, 
      WorkingPrecision -&gt; prec] - 
     RandomVariate[ExponentialDistribution[b], n, 
      WorkingPrecision -&gt; prec] + 
     RandomVariate[NormalDistribution[m, s], n, 
      WorkingPrecision -&gt; prec];
  
  (* Fitting: This uses Mean, central moments 2 and 3 and 4th cumulant \
  but it often does not provide a solution *)
  
  nDistParam[data_] := Module[{mn, vv, m3, k4, al, be, m, si},
        mn = Mean[data];
        vv = CentralMoment[data, 2];
        m3 = CentralMoment[data, 3];
        k4 = Cumulant[data, 4];
        al = 
      ConditionalExpression[
       Root[864 - 864 m3 #1^3 - 216 k4 #1^4 + 648 m3^2 #1^6 + 
          36 k4^2 #1^8 - 216 m3^3 #1^9 + (-2 k4^3 + 27 m3^4) #1^12 &amp;, 
        2], k4 &gt; Root[-27 m3^4 + 4 #1^3 &amp;, 1]];
        be = ConditionalExpression[
  
       Root[2 Root[
             864 - 864 m3 #1^3 - 216 k4 #1^4 + 648 m3^2 #1^6 + 
               36 k4^2 #1^8 - 
               216 m3^3 #1^9 + (-2 k4^3 + 27 m3^4) #1^12 &amp;, 
             2]^3 + (-2 + 
             m3 Root[
                864 - 864 m3 #1^3 - 216 k4 #1^4 + 648 m3^2 #1^6 + 
                  36 k4^2 #1^8 - 
                  216 m3^3 #1^9 + (-2 k4^3 + 27 m3^4) #1^12 &amp;, 
                2]^3) #1^3 &amp;, 1], k4 &gt; Root[-27 m3^4 + 4 #1^3 &amp;, 1]];
        m = mn - 1/al + 1/be;
        si = 
      Sqrt[Abs[-al^-2 - be^-2 + vv ]];(*Ensure positive*)
        {al, 
      be, m, si}];
  
  nDistLL = 
    Compile[{a, b, m, s, {x, _Real, 1}}, 
     Total[Log[
       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)])]](*, CompilationTarget-&gt;"C", 
     RuntimeAttributes-&gt;{Listable}, Parallelization-&gt;True*)];
  
  nlloglike[data_, a_?NumericQ, b_?NumericQ, m_?NumericQ, s_?NumericQ] := 
    nDistLL[a, b, m, s, data];
  
  nFit[data_] := Module[{a, b, m, s, a0, b0, m0, s0, res},
  
        (* So far have not found a good way to quickly estimate a and \
  b.  Starting assumption is that they both = 2,then m0 ~= 
     Mean and s0 ~= 
     StandardDeviation it seems to work better if a and b are not the \
  same at start. *)
  
     {a0, b0, m0, s0} = nDistParam[data];(*may give Undefined values*)
  
       If[! (VectorQ[{a0, b0, m0, s0}, NumericQ] &amp;&amp; 
         VectorQ[{a0, b0, s0}, # &gt; 0 &amp;]),
              m0 = Mean[data];
              s0 = StandardDeviation[data];
              a0 = 1;
              b0 = 2;];
     res = {a, b, m, s} /. 
       FindMaximum[
         nlloglike[data, Abs[a], Abs[b], m,  
          Abs[s]], {{a, a0}, {b, b0}, {m, m0}, {s, s0}},
                 Method -&gt; "PrincipalAxis"][[2]];
        {Abs[res[[1]]], Abs[res[[2]]], res[[3]], Abs[res[[4]]]}];
  
  nFit[data_, {a0_, b0_, m0_, s0_}] := Module[{a, b, m, s, res},
        res = {a, b, m, s} /. 
       FindMaximum[
         nlloglike[data, Abs[a], Abs[b], m, 
          Abs[s]], {{a, a0}, {b, b0}, {m, m0}, {s, s0}},
                 Method -&gt; "PrincipalAxis"][[2]];
        {Abs[res[[1]]], Abs[res[[2]]], res[[3]], Abs[res[[4]]]}];
  
  dDist /: PDF[dDist[a_, b_, m_, s_], x_] := 
    PDF[nDist[a, b, m, s], Log[x]]/x;
  dDist /: CDF[dDist[a_, b_, m_, s_], x_] := 
    CDF[nDist[a, b, m, s], Log[x]];
  dDist /: EstimatedDistribution[data_, dDist[a_, b_, m_, s_]] := 
    dDist[Sequence @@ nFit[Log[data]]];
  dDist /: EstimatedDistribution[data_, 
     dDist[a_, b_, m_, 
      s_], {{a_, a0_}, {b_, b0_}, {m_, m0_}, {s_, s0_}}] := 
    dDist[Sequence @@ nFit[Log[data], {a0, b0, m0, s0}]];
  dDist /: Quantile[dDist[a_, b_, m_, s_], p_] := 
   Module[{x}, x /. FindRoot[CDF[dDist[a, b, m, s], x] == p, {x, s}]] /;
     0 &lt; p &lt; 1
  dDist /: Quantile[dDist[a_, b_, m_, s_], p_] :=  
   Module[{x}, 
     x /. FindRoot[ CDF[dDist[a, b, m, s], x] == #, {x, s}] &amp; /@ p] /; 
    VectorQ[p, 0 &lt; # &lt; 1 &amp;]
  dDist /: Quantile[dDist[a_, b_, m_, s_], p_] := -Infinity /; p == 0
  dDist /: Quantile[dDist[a_, b_, m_, s_], p_] := Infinity /; p == 1
  dDist /: DistributionDomain[dDist[a_, b_, m_, s_]] := 
   Interval[{0, Infinity}]
  dDist /: DistributionParameterQ[dDist[a_, b_, m_, s_]] := ! 
    TrueQ[Not[Element[{a, b, s, m}, Reals] &amp;&amp; a &gt; 0 &amp;&amp; b &gt; 0 &amp;&amp; s &gt; 0]]
  dDist /: DistributionParameterAssumptions[dDist[a_, b_, m_, s_]] := 
   Element[{a, b, s, m}, Reals] &amp;&amp; a &gt; 0 &amp;&amp; b &gt; 0 &amp;&amp; s &gt; 0
  dDist /: Random`DistributionVector[dDist[a_, b_, m_, s_], n_, prec_] :=
     Exp[RandomVariate[ExponentialDistribution[a], n, 
       WorkingPrecision -&gt; prec] - 
         RandomVariate[ExponentialDistribution[b], n, 
       WorkingPrecision -&gt; prec] + 
      RandomVariate[NormalDistribution[m, s], n, 
       WorkingPrecision -&gt; prec]];
  </code></pre>
  
  <p>This enables me to fit distribution parameters and generate <strong>PDF's</strong> and <strong>CDF's</strong>.  An example of the plots:</p>
  
  <pre><code>Plot[PDF[dDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
   PlotRange -&gt; All]
  Plot[CDF[dDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
   PlotRange -&gt; All]
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/0ynOS.png" alt="enter image description here"></p>
  
  <p>Now I've defined a <code>function</code> to calculate mean residual life (see <a href="https://math.stackexchange.com/questions/104494/mean-residual-life-clarification-needed">this question</a> for an explanation).</p>
  
  <pre><code>MeanResidualLife[start_, dist_] := 
   NExpectation[X \[Conditioned] X &gt; start, X \[Distributed] dist] - 
    start
  MeanResidualLife[start_, limit_, dist_] := 
   NExpectation[X \[Conditioned] start &lt;= X &lt;= limit, 
     X \[Distributed] dist] - start
  </code></pre>
  
  <p>The first of these that doesn't set a limit as in the second takes a long time to calculate, but they both work.</p>
  
  <p>Now I need to find the minimum of the <code>MeanResidualLife</code> function for the same distribution (or some variation of it) or minimize it.</p>
  
  <p>I've tried a number of variations on this:</p>
  
  <pre><code>FindMinimum[MeanResidualLife[x, dDist[3.77, 1.34, -2.65, 0.40]], x]
  FindMinimum[MeanResidualLife[x, 1, dDist[3.77, 1.34, -2.65, 0.40]], x]
  
  NMinimize[{MeanResidualLife[x, dDist[3.77, 1.34, -2.65, 0.40]], 
    0 &lt;= x &lt;= 1}, x]
  NMinimize[{MeanResidualLife[x, 1, dDist[3.77, 1.34, -2.65, 0.40]], 0 &lt;= x &lt;= 1}, x]
  </code></pre>
  
  <p>These either seem to run forever or run into: </p>
  
  <blockquote>
    <p>Power::infy : Infinite expression 1/ 0. encountered. >></p>
  </blockquote>
  
  <p>The <code>MeanResidualLife</code> function applied to a simpler but similarly shaped distribution shows that it has a single minimum:</p>
  
  <pre><code>Plot[PDF[LogNormalDistribution[1.75, 0.65], x], {x, 0, 30}, 
   PlotRange -&gt; All]
  Plot[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], {x, 0, 
    30},
   PlotRange -&gt; {{0, 30}, {4.5, 8}}]
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/4IpRL.png" alt="enter image description here"></p>
  
  <p>Also both:</p>
  
  <pre><code>FindMinimum[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], x]
  FindMinimum[MeanResidualLife[x, 30, LogNormalDistribution[1.75, 0.65]], x]
  </code></pre>
  
  <p>give me answers (if with a bunch of messages first) when used with the <code>LogNormalDistribution</code>.</p>
  
  <p>Any thoughts on how to get this to work for the custom distribution described above?</p>
  
  <p>Do I need to add constraints or options?</p>
  
  <p>Do I need to define something else in the definitions of the custom distributions?</p>
  
  <p>Maybe the <code>FindMinimum</code> or <code>NMinimize</code> just need to run longer (I've run them nearly an hour to no avail).  If so do I just need some way to speed up finding the minimum of the function?  Any suggestions on how?</p>
  
  <p>Does <code>Mathematica</code> have another way to do this?</p>
  
  <p><strong>Added 9 Feb 5:50PM EST:</strong> </p>
  
  <p>Anyone can download <strong>Oleksandr Pavlyk's</strong> presentation about creating distributions in Mathematica from the Wolfram Technology Conference 2011 workshop 'Create Your Own Distribution' <a href="http://www.wolfram.com/events/technology-conference/2011/mathematics-and-statistics.html" rel="nofollow noreferrer">here</a>.  The downloads include the notebook, <code>'ExampleOfParametricDistribution.nb'</code> that seems to lays out all the pieces required to create a distribution that one can use like the distributions that come with Mathematica.</p>
  
  <p>It may supply some of the answer.</p>