Вопрос по probability, plot, calculus, wolfram-mathematica, mathematica-8 – Минимизация NExpectation для пользовательского дистрибутива в Mathematica

237

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

Расчет ожидания для пользовательского распределения в 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.

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

Не эксперт по Mathematica, но я сталкивался с подобными проблемами в других местах. Кажется, у вас возникают проблемы, когда ваш домен начинается с 0. Попробуйте начать с 0.1 и выше и посмотрите, что произойдет. Makketronix
@Makketronix - Спасибо за это. Забавная синхронность, учитывая, что я начал возвращаться к этому через 3 года. Jagra
Я не уверен, что смогу помочь вам, но вы могли бы попытаться спросить уMathematica-специфичный стекопоток, Удачи! Olivia Stork
Есть куча статей об этом наzbmath.org Поиск ожиданий Ivan V

Ваш Ответ

1   ответ
11

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

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

Сначала давайте посмотрим, чтоMeanResidualLife есть, как вы это определили.NExpectation или жеExpectation вычислитьожидаемое значение, Для ожидаемого значения нам нужно только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, я думаю, нам нужно нормализовать результат делением на интеграл PDF в этом диапазоне. Таким образом, мое предположение о левом ожидаемом значении

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

И для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, {start -> 0.0584312}}

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

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

+1 - только что видел, так что мне нужно будет проработать это, но я думаю, что способ, которым вы разделили проблему на решаемые шаги, имеет большой смысл. Кроме того, сюжет вашей функции MRL, безусловно, выглядит на месте. Большое спасибо, я вернусь к этому, как только смогу найти время для изучения вашего ответа. Jagra

Похожие вопросы