Минимизация NExpectation для пользовательского распределения в Mathematica
Это относится к более раннему вопросу еще в июне:
Вычисление ожидания для пользовательского распределения в 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]
![enter image description here]()
Теперь я определил a 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}}]
![enter image description here]()
Также оба:
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:
Любой может скачать презентацию Oleksandr Pavlyk's о создании дистрибутивов в Mathematica с мастерской Wolfram Technology Conference 2011 "Создайте собственное расписание" здесь. Загрузки включают в себя блокнот, 'ExampleOfParametricDistribution.nb', который, как представляется, содержит все части, необходимые для создания дистрибутива, который можно использовать, например, в дистрибутивах, поставляемых с Mathematica.
Он может предоставить часть ответа.
Ответы
Ответ 1
Насколько я вижу, проблема (как вы уже писали), что 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
Plot[pdf2[3.77, 1.34, -2.65, 0.40, x], {x, 0, .3}]
![Plot of PDF]()
Теперь к ожидаемому значению. Если я правильно ее понимаю, мы должны интегрировать 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]
![Plot of x * PDF]()
и ожидаемое значение
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}]
![Plot of Mean Residual Life]()
Выглядит правдоподобно, но я не эксперт. Итак, мы хотим минимизировать его, т.е. Найти start, для которого эта функция является локальным минимумом. Кажется, что минимум составляет около 0,05, но можно найти более точное значение, начиная с этой догадки
FindMinimum[MRL[start], {start, 0.05}]
и после некоторых ошибок (ваша функция не определена ниже 0, поэтому я думаю, что минимизатор немного замалчивает в этой запретной области), мы получаем
{0.0418137, {start → 0.0584312}}
Таким образом, оптимальный должен быть на start = 0.0584312 со средним остаточным сроком 0.0418137.
Я не знаю, правильно ли это, но кажется правдоподобным.
