在Mathematica中最小化NExpectation以获得自定义分布
这涉及到6月份的一个早期问题：
计算Mathematica中自定义分布的期望
我有一个自定义的混合分布定义使用第二个自定义分布按照@Sasha在过去一年的答案中讨论的@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有另一种方法来做到这一点吗？
2月9日下午5:50加入EST：
任何人都可以在这里下载Wolfram Technology Conference 2011研讨会的“创建您自己的分销”研讨会，下载Oleksandr Pavlyk关于在Mathematica中创建分销的演示。 下载内容包括笔记本'ExampleOfParametricDistribution.nb' ，它似乎列出了创建一个可以像Mathematica附带的发行版一样使用的发行版所需的所有部分。
它可能会提供一些答案。
就我所见，问题是（正如你已经写过的那样），即使对于单个评估， 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 。
我不知道这是否正确，但似乎是合理的。
链接地址: http://www.djcxy.com/p/6055.html上一篇: Minimizing NExpectation for a custom distribution in Mathematica