Minimalizowanie NExpectation dla niestandardowej dystrybucji w Mathematica

Odnosi się to do wcześniejszego pytania z czerwca:

Obliczanie oczekiwanego rozkładu niestandardowego w Mathematica

Mam niestandardową dystrybucję mieszaną zdefiniowaną za pomocą drugiej dystrybucji niestandardowej zgodnie z wytycznymi omówionymi @Sashaw wielu odpowiedziach w ciągu ostatniego roku.

Kod definiujący dystrybucje wygląda następująco:

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]];

To pozwala mi dopasować parametry dystrybucji i wygenerować pliki PDF i CDF . Przykład fabuły:

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]

Teraz zdefiniowałem a functiondo obliczenia średniego pozostałego życia (zobacz wyjaśnienie w tym pytaniu ).

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

Pierwszy z nich, który nie określa limitu, jak w drugim, zajmuje dużo czasu, ale oba działają.

Teraz muszę znaleźć minimum MeanResidualLife funkcji dla tego samego rozkładu (lub jakiejś jego odmiany) lub zminimalizować.

Wypróbowałem wiele odmian tego:

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]

Te wydają się działać wiecznie lub napotkać:

Power :: infy: Napotkano nieskończone wyrażenie 1 / 0. >>

MeanResidualLifeFunkcja stosowane do prostych, ale podobnie ukształtowanych pokazuje rozkład, że ma co najmniej jeden:

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}}]

Również oba:

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

daj mi odpowiedzi (jeśli najpierw jest kilka wiadomości), jeśli używasz go razem z LogNormalDistribution.

Czy są jakieś przemyślenia na temat tego, jak to zrobić w przypadku niestandardowej dystrybucji opisanej powyżej?

Czy muszę dodawać ograniczenia lub opcje?

Czy muszę definiować coś jeszcze w definicjach niestandardowych dystrybucji?

Może FindMinimumlubNMinimize po prostu muszę biec dłużej (prowadzę je prawie godzinę bezskutecznie). Jeśli tak, to czy potrzebuję jakiegoś sposobu, aby przyspieszyć znalezienie minimum funkcji? Wszelkie sugestie jak?

Czy Mathematicama na to inny sposób?

Dodano 9 lutego 17:50 EST:

Każdy może pobrać prezentację Oleksandra Pavlyka na temat tworzenia dystrybucji w Mathematica z warsztatów Wolfram Technology Conference 2011 „Stwórz własną dystrybucję” tutaj . Pliki do pobrania obejmują notatnik, 'ExampleOfParametricDistribution.nb'który wydaje się określać wszystkie elementy wymagane do utworzenia dystrybucji, z której można korzystać podobnie jak dystrybucje dostarczane z Mathematica.

Może dostarczyć niektórych odpowiedzi.

O ile mi wiadomo, problemem jest (jak już pisałeś), którego MeanResidualLifeobliczenie zajmuje dużo czasu, nawet w przypadku pojedynczej oceny. TerazFindMinimum lub podobne funkcje próbują znaleźć minimum dla funkcji. Znalezienie minimum wymaga ustawienia pierwszej pochodnej funkcji zero i rozwiązania rozwiązania. Ponieważ twoja funkcja jest dość skomplikowana (i prawdopodobnie nie można jej rozróżnić), drugą możliwością jest minimalizacja numeryczna, która wymaga wielu ocen twojej funkcji. Ergo, to jest bardzo, bardzo wolne.

Sugerowałbym wypróbowanie go bez magii Mathematica.

Najpierw zobaczmy, co to MeanResidualLifejest, jak to zdefiniowałeś. NExpectationlub Expectationobliczyć oczekiwaną wartość . Aby uzyskać oczekiwaną wartość, potrzebujemy tylko PDFtwojej dystrybucji. Wyodrębnijmy go z powyższej definicji na proste funkcje:

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;

Jeśli wydrukujemy pdf2, będzie wyglądać dokładnie tak samo jak twoja fabuła

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

Teraz do oczekiwanej wartości. Jeśli dobrze to rozumiem, musimy zintegrować x * pdf[x]od -infdo +infdla normalnej oczekiwanej wartości.

x * pdf[x] wygląda jak

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

a oczekiwana wartość to

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

Ale ponieważ chcesz oczekiwaną wartość między a starti +infmusimy zintegrować w tym zakresie, a ponieważ PDF nie integruje się już z 1 w tym mniejszym przedziale, myślę, że musimy znormalizować wynik dzieląc przez całkę z pliku PDF w ten zakres. Tak więc przypuszczam, że wartość oczekiwana po lewej stronie jest

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]}]

A dla tego, MeanResidualLifeco startod niego odejmujesz , dając

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

Które działki jak

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

Wygląda na wiarygodne, ale nie jestem ekspertem. W końcu chcemy to zminimalizować, tzn. Znaleźć, startdla którego ta funkcja jest lokalnym minimum. Minimum wydaje się wynosić około 0,05, ale znajdźmy dokładniejszą wartość, zaczynając od tego przypuszczenia

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

i po kilku błędach (twoja funkcja nie jest zdefiniowana poniżej 0, więc myślę, że minimalizator trochę zaczepia w tym zabronionym regionie) otrzymujemy

{0,0418137, {start -> 0,0584312}}

Zatem optymalne powinno być przy start = 0.0584312średnim okresie resztkowym wynoszącym 0.0418137.

Nie wiem, czy to prawda, ale wydaje się prawdopodobne.