最大限度地减lessMathematica中的自定义分配的NExpectation

这涉及到六月份以前的一个问题：

在Mathematica中计算自定义分布的期望

我有一个自定义混合分布定义使用第二个自定义分布按照@Sasha在过去一年的答案中讨论的@Sasha行。

定义分布的代码如下： 

 nDist /: CharacteristicFunction[nDist[a_, b_, m_, s_], t_] := (ab E^(I mt - (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)) ab (E^(a (m + (as^2)/2 - x)) Erfc[(m + as^2 - x)/(Sqrt[2] s)] + E^(b (-m + (bs^2)/2 + x)) Erfc[(-m + bs^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

第一个没有像第二个那样设定限制需要很长时间来计算，但是它们都起作用。

现在我需要find相同分布（或其一些变化）的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：遇到无限expression1 / 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只需要运行更长的时间（我已经运行了近一个小时无济于事）。 如果是这样，我只需要一些方法来加快find函数的最小值？ 任何build议如何？

Mathematica有另一种方法来做到这一点？

2月9日5:50 PM EST：

任何人都可以从这里下载Wolfram Technology Conference 2011研讨会“创build你自己的发行版”，下载Oleksandr Pavlyk关于在Mathematica中创build发行版的演示。 下载内容包括笔记本'ExampleOfParametricDistribution.nb' ，它似乎列出了创build一个可以像Mathematica附带的发行版一样使用的所有发行版。

它可能会提供一些答案。

就我所见，问题是（正如你已经写过的那样），即使是单一评估， MeanResidualLife需要很长的时间来计算。 现在， FindMinimum或类似的函数试图find函数的最小值。 find最小值需要设置函数零点的一阶导数并求解。 由于你的函数相当复杂（可能不可微分），第二种可能性是做一个数值最小化，这需要你的函数的许多评价。 埃尔戈，这是非常非常缓慢的。

我会build议尝试没有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}]

PDF的地图

现在到预期的价值。 如果我理解正确，我们必须从-inf到+inf整合x * pdf[x]以获得正常的期望值。

x * pdf[x]看起来像 

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

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

平均剩余寿命图

看起来似乎合理，但我不是专家。 所以最后我们想把它最小化，即find这个函数是一个局部最小值的开始。 最小值似乎在0.05左右，但是从这个猜测开始，我们可以find一个更精确的值 

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

和一些错误之后（你的函数没有被定义在0以下，所以我想最小化者在那个被禁止的区域内稍微摸了一下），我们得到

{0.0418137，{start – > 0.0584312}}

所以最佳值应该是在start = 0.0584312 ，平均剩余寿命为0.0418137 。

我不知道这是否正确，但似乎是合理的。

Interesting Posts