Это относится к более раннему вопросу еще в июне:

Расчет ожидания для пользовательского распределения в Mathematica

У меня есть пользовательский смешанный дистрибутив, определенный с помощью второго пользовательского дистрибутива в соответствии с тем, как@Sasha в ряде ответов за прошедший год.

Код, определяющий распределение:

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

Это позволяет мне соответствовать параметрам распределения и генерироватьPDF-х а такжеCDF-х, Пример участков:

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]

Теперь я определилfunction рассчитать средний остаточный срок службы (см.этот вопрос для объяснения).

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

Первый из них, который не устанавливает ограничения, как во втором, занимает много времени для расчета, но они оба работают.

Теперь мне нужно найти минимумMeanResidualLife функция для того же самого распределения (или некоторого изменения этого) или минимизируйте это.

Я пробовал несколько вариантов этого:

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]

Кажется, что они либо навсегда, либо сталкиваются с:

Power :: infy: встречается бесконечное выражение 1 / 0. >>

MeanResidualLife Функция, примененная к более простому, но похожему по форме распределению, показывает, что она имеет единственный минимум:

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

Также оба:

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

дайте мне ответы (если с кучей сообщений в первую очередь) при использовании сLogNormalDistribution.

Любые мысли о том, как заставить это работать для пользовательского дистрибутива, описанного выше?

Нужно ли добавлять ограничения или параметры?

Нужно ли мне определять что-то еще в определениях пользовательских дистрибутивов?

Может бытьFindMinimum или жеNMinimize просто нужно бежать дольше (я бегал их почти час безрезультатно). Если так, то мне просто нужен какой-то способ ускорить поиск минимума функции? Любые предложения о том, как?

Есть лиMathematica есть другой способ сделать это?

Добавлено 9 фев. 17:50 EST:

Любой может скачатьАлександра Павлика презентация о создании дистрибутивов в Mathematica на семинаре Wolfram Technology Conference 2011 «Создайте свой дистрибутив»Вот, Загрузки включают ноутбук,'ExampleOfParametricDistribution.nb' Кажется, в нем изложены все элементы, необходимые для создания дистрибутива, который можно использовать, например, дистрибутивы, поставляемые с Mathematica.

Это может дать некоторые ответы.