Minimizando o NExpectation para uma distribuição personalizada no Mathematica

Isso se refere a uma pergunta anterior de junho:

Calcular a expectativa de uma distribuição personalizada no Mathematica

Tenho uma distribuição mista personalizada definida usando uma segunda distribuição personalizada seguindo as linhas discutidas por@Sasha em várias respostas no ano passad

Code definindo as distribuições a seguir:

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}] & /@ p] /; 
  VectorQ[p, 0 < # < 1 &]
nDist /: Quantile[nDist[a_, b_, m_, s_], p_] := 
 Module[{x}, x /. FindRoot[CDF[nDist[a, b, m, s], x] == p, {x, m}]] /;
   0 < p < 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] && a > 0 && b > 0 && s > 0]]
nDist /: DistributionParameterAssumptions[nDist[a_, b_, m_, s_]] := 
 Element[{a, b, s, m}, Reals] && a > 0 && b > 0 && s > 0
nDist /: Random`DistributionVector[nDist[a_, b_, m_, s_], n_, prec_] :=

    RandomVariate[ExponentialDistribution[a], n, 
    WorkingPrecision -> prec] - 
   RandomVariate[ExponentialDistribution[b], n, 
    WorkingPrecision -> prec] + 
   RandomVariate[NormalDistribution[m, s], n, 
    WorkingPrecision -> 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 &, 
      2], k4 > Root[-27 m3^4 + 4 #1^3 &, 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 &, 
           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 &, 
              2]^3) #1^3 &, 1], k4 > Root[-27 m3^4 + 4 #1^3 &, 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->"C", 
   RuntimeAttributes->{Listable}, Parallelization->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] && 
       VectorQ[{a0, b0, s0}, # > 0 &]),
            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 -> "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 -> "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 < p < 1
dDist /: Quantile[dDist[a_, b_, m_, s_], p_] :=  
 Module[{x}, 
   x /. FindRoot[ CDF[dDist[a, b, m, s], x] == #, {x, s}] & /@ p] /; 
  VectorQ[p, 0 < # < 1 &]
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] && a > 0 && b > 0 && s > 0]]
dDist /: DistributionParameterAssumptions[dDist[a_, b_, m_, s_]] := 
 Element[{a, b, s, m}, Reals] && a > 0 && b > 0 && s > 0
dDist /: Random`DistributionVector[dDist[a_, b_, m_, s_], n_, prec_] :=
   Exp[RandomVariate[ExponentialDistribution[a], n, 
     WorkingPrecision -> prec] - 
       RandomVariate[ExponentialDistribution[b], n, 
     WorkingPrecision -> prec] + 
    RandomVariate[NormalDistribution[m, s], n, 
     WorkingPrecision -> prec]];

Isto me permite ajustar parâmetros de distribuição e gerar PDF's e CDF's. Um exemplo das parcelas:

Plot[PDF[dDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
 PlotRange -> All]
Plot[CDF[dDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
 PlotRange -> All]

Agora eu defini umfunction para calcular a vida útil média residual (consulteessa questã para uma explicação

MeanResidualLife[start_, dist_] := 
 NExpectation[X \[Conditioned] X > start, X \[Distributed] dist] - 
  start
MeanResidualLife[start_, limit_, dist_] := 
 NExpectation[X \[Conditioned] start <= X <= limit, 
   X \[Distributed] dist] - start

O primeiro desses que não define um limite, pois no segundo leva muito tempo para calcular, mas ambos funciona

Agora preciso encontrar o mínimo deMeanResidualLifeunção para a mesma distribuição (ou alguma variação dela) ou minimizá-l

Eu tentei várias variações sobre isso:

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 <= x <= 1}, x]
NMinimize[{MeanResidualLife[x, 1, dDist[3.77, 1.34, -2.65, 0.40]], 0 <= x <= 1}, x]

Estes parecem durar para sempre ou deparar-se com:

Power :: infy: expressão infinita 1 / 0. encontrada. >>

OMeanResidualLife função @ aplicada a uma distribuição mais simples, mas com formato semelhante, mostra que ela possui um único mínimo:

Plot[PDF[LogNormalDistribution[1.75, 0.65], x], {x, 0, 30}, 
 PlotRange -> All]
Plot[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], {x, 0, 
  30},
 PlotRange -> {{0, 30}, {4.5, 8}}]

Também ambos:

FindMinimum[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], x]
FindMinimum[MeanResidualLife[x, 30, LogNormalDistribution[1.75, 0.65]], x]

me dê respostas (se houver várias mensagens primeiro) quando usado com oLogNormalDistribution.

Alguma idéia de como fazer isso funcionar para a distribuição personalizada descrita acim

Preciso adicionar restrições ou opções?

reciso definir outra coisa nas definições das distribuições personalizada

Talvez oFindMinimum ouNMinimize só preciso correr mais tempo (já corri quase uma hora sem sucesso). Em caso afirmativo, só preciso de uma maneira de acelerar a localização do mínimo da função? Alguma sugestão de como?

FazMathematica tem outra maneira de fazer isso?

Added 9 de fev 17:50 EST:

ualquer pessoa pode fazer o downloadOleksandr Pavlyk's apresentação sobre a criação de distribuições no Mathematica a partir do workshop Wolfram Technology Conference 2011 'Crie sua própria distribuição'aqu. Os downloads incluem o notebook,'ExampleOfParametricDistribution.nb' que parece apresentar todas as peças necessárias para criar uma distribuição que se possa usar como as distribuições que acompanham o Mathematic

Pode fornecer algumas das resposta

questionAnswers(2)

yourAnswerToTheQuestion