Это относится к более раннему вопросу еще в июне:
Расчет ожидания для пользовательского распределения в 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.
Это может дать некоторые ответы.
Ответы:
Насколько я понимаю, проблема в том, что (как вы уже писали),
MeanResidualLife
на вычисление уходит много времени, даже для одной оценки. ТеперьFindMinimum
аналогичные функции пытаются найти минимум функции. Для нахождения минимума требуется либо установить первую производную функции ноль, либо найти решение. Поскольку ваша функция довольно сложна (и, вероятно, не дифференцируема), вторая возможность состоит в том, чтобы выполнить числовую минимизацию, которая требует много оценок вашей функции. Поэтому очень и очень медленно.Я бы предложил попробовать это без магии Mathematica.
Сначала давайте посмотрим, что это
MeanResidualLife
такое, как вы это определили.NExpectation
илиExpectation
рассчитать ожидаемое значение . Для ожидаемого значения нам нужен толькоPDF
ваш дистрибутив. Давайте выделим это из вашего определения выше в простые функции:Если мы строим pdf2, он выглядит точно так же, как ваш сюжет
Теперь к ожидаемому значению. Если я правильно понимаю, мы должны интегрировать
x * pdf[x]
от-inf
к+inf
для нормального ожидаемого значения.x * pdf[x]
выглядит каки ожидаемое значение
Но так как вы хотите , ожидаемое значение между
start
и+inf
нам нужно было интегрировать в этом диапазоне, и, поскольку PDF больше не интегрируется в 1 в этом меньшем интервале, я думаю, что мы должны нормализовать результат делением на интеграл PDF в этот диапазон. Таким образом, мое предположение о левом ожидаемом значенииИ для этого
MeanResidualLife
вы вычитаетеstart
, даваяКакие участки как
Выглядит правдоподобно, но я не эксперт. Итак, наконец, мы хотим минимизировать его, т.е. найти,
start
для которого эта функция является локальным минимумом. Минимум кажется около 0,05, но давайте найдем более точное значение, исходя из этого предположенияи после некоторых ошибок (ваша функция не определена ниже 0, так что я думаю, минимизатор немного ткнет в эту запрещенную область), мы получаем
{0.0418137, {start -> 0.0584312}}
Таким образом, оптимум должен быть при
start = 0.0584312
средней остаточной жизни0.0418137
.Я не знаю, правильно ли это, но это кажется правдоподобным.
источник