Минимизация NExpectation для распределения пользовательских настроек в системе Mathematica

Question

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

Calculating expectation for a custom distribution in 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 для расчета среднего остаточного срока службы (см this question для объяснения).

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] 

Они либо как запустить навсегда или запустить в:

Мощность :: инфы: Infinite выражения 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 февраля 5:50 вечера EST:

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

Он может предоставить некоторые ответы.

Answer 1

Насколько я вижу, проблема (как вы уже писали), что MeanResidualLife занимает много времени, чтобы вычислить, даже для одной оценки. Теперь FindMinimum или аналогичные функции пытаются найти минимум функции. Для нахождения минимума требуется либо установить первую производную функции нуль, либо решить ее для решения. Поскольку ваша функция довольно сложная (и, вероятно, недифференцируемая), вторая возможность заключается в численной минимизации, что требует многих оценок вашей функции. Эрго, это очень медленно.

Я предлагаю попробовать его без магии Mathematica.

Сначала посмотрим, что такое MeanResidualLife, как вы его определили. NExpectation или Expectation вычислить expected value. Для ожидаемого значения нам нужен только PDF вашего дистрибутива. Давайте извлечь его из определения выше в простые функции:

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; 

Если мы наносим PDF2 он выглядит точно так же, как ваш участок

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

Теперь к ожидаемому значению. Если я правильно понимаю, мы должны интегрировать x * pdf[x] с -inf в +inf для нормального ожидаемого значения.

x * pdf[x] выглядит

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

и ожидаемое значение

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

Но так как вы хотите, ожидаемое значение между start и +inf нам нужно интегрировать в этом диапазоне , и так как PDF тогда больше не интегрируется в 1 в этом меньшем интервале, я думаю, мы hav e для нормализации результата делятся на интеграл PDF в этом диапазоне.Так что я думаю для левой переплете ожидаемого значения является

А для MeanResidualLife вычесть start из него, давая

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

, который участки, как

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

Выглядит правдоподобно, но я не эксперт. Поэтому, наконец, мы хотим минимизировать его, т. Е. Найти start, для которого эта функция является локальным минимумом. Минимум, кажется, около 0,05, но давайте найдем более точное значение, начинающееся с этой догадкой

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

и после некоторых ошибок (ваша функция не определена ниже 0, так что я думаю, минимизирующие тыкаю немного в том, что запрещено область) мы получаем

{0,0418137, {старт -> 0,0584312}}

Таким образом, оптимальные должен быть start = 0.0584312 со средней остаточной жизнью 0.0418137.

Я не знаю, правильно ли это, но кажется правдоподобным.