为Mathematica中的自定义分布最小化NExpectation [英] Minimizing NExpectation for a custom distribution in Mathematica
问题描述
这与6月份的一个早期问题有关:
This relates to an earlier question from back in June:
在过去的一年中，我根据第二个自定义分布定义了一个自定义混合分布，该分布遵循@Sasha在很多答案中讨论的内容.
I have a custom mixed distribution defined using a second custom distribution following along the lines discussed by @Sasha in a number of answers over the past year.
定义分布的代码如下:
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的.情节的例子:
This enables me to fit distribution parameters and generate PDF's and CDF's. An example of the plots:
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来计算平均剩余寿命(请参阅此问题进行解释).
Now I've defined a function to calculate mean residual life (see this question for an explanation).
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
其中第一个没有设置限制，而第二个设置花费很长时间，但是它们都可以工作.
The first of these that doesn't set a limit as in the second takes a long time to calculate, but they both work.
现在，我需要找到相同分布(或它的某些变化)的MeanResidualLife函数的最小值或将其最小化.
Now I need to find the minimum of the MeanResidualLife function for the same distribution (or some variation of it) or minimize it.
我对此进行了多种尝试:
I've tried a number of variations on this:
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]
这些似乎永远消失或碰到:
These either seem to run forever or run into:
Power :: infy:遇到无限表达式1/0. >>
Power::infy : Infinite expression 1/ 0. encountered. >>
应用于更简单但形状相似的分布的MeanResidualLife函数表明它只有一个最小值:
The MeanResidualLife function applied to a simpler but similarly shaped distribution shows that it has a single minimum:
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一起使用时，给我答案(如果首先显示一堆消息).
give me answers (if with a bunch of messages first) when used with the LogNormalDistribution.
是否有任何想法如何使它适用于上述自定义发行版?
Any thoughts on how to get this to work for the custom distribution described above?
我需要添加约束或选项吗?
Do I need to add constraints or options?
我需要在自定义发行版的定义中定义其他内容吗?
Do I need to define something else in the definitions of the custom distributions?
也许FindMinimum或NMinimize只需要运行更长的时间(我已经将它们运行了将近一个小时而无济于事).如果是这样，我是否仅需要某种方法来加快查找功能的最小值?有什么建议吗?
Maybe the FindMinimum or NMinimize just need to run longer (I've run them nearly an hour to no avail). If so do I just need some way to speed up finding the minimum of the function? Any suggestions on how?
Mathematica是否有另一种方法?
于2月9日美国东部时间下午5:50添加:
任何人都可以从Wolfram技术大会2011研讨会创建自己的分布"中下载 Oleksandr Pavlyk关于在Mathematica中创建分布的演示文稿
Anyone can download Oleksandr Pavlyk's presentation about creating distributions in Mathematica from the Wolfram Technology Conference 2011 workshop 'Create Your Own Distribution' here. The downloads include the notebook, 'ExampleOfParametricDistribution.nb' that seems to lays out all the pieces required to create a distribution that one can use like the distributions that come with Mathematica.
它可能会提供一些答案.
It may supply some of the answer.
推荐答案
据我所知，问题是(正如您已经写过的)，即使对于单个评估，MeanResidualLife也需要花费很长的时间来计算.现在，FindMinimum或类似的函数尝试查找该函数的最小值.找到最小值需要将函数的一阶导数设置为零并求解.由于您的函数非常复杂(并且可能不可微)，因此第二种可能性是进行数值最小化，这需要对函数进行多次评估.嗯，这非常非常慢.
As far as I see, the problem is (as you already wrote), that MeanResidualLife takes a long time to compute, even for a single evaluation. Now, the FindMinimum or similar functions try to find a minimum to the function. Finding a minimum requires either to set the first derivative of the function zero and solve for a solution. Since your function is quite complicated (and probably not differentiable), the second possibility is to do a numerical minimization, which requires many evaluations of your function. Ergo, it is very very slow.
我建议您在没有Mathematica魔术的情况下尝试使用它.
I'd suggest to try it without Mathematica magic.
首先，让我们看看定义的MeanResidualLife是什么. NExpectation或Expectation计算期望值.
对于期望值，我们只需要您的发行版的PDF.让我们从上面的定义中将其提取为简单的函数:
First let's see what the MeanResidualLife is, as you defined it. NExpectation or Expectation compute the expected value.
For the expected value, we only need the PDF of your distribution. Let's extract it from your definition above into simple functions:
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，则其外观与您的绘制完全相同
If we plot pdf2 it looks exactly as your Plot
Plot[pdf2[3.77, 1.34, -2.65, 0.40, x], {x, 0, .3}]
现在达到期望值.如果我理解正确，我们必须将x * pdf[x]从-inf集成到+inf以获得正常的期望值.
Now to the expected value. If I understand it correctly we have to integrate x * pdf[x] from -inf to +inf for a normal expected value.
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的整数.所以我对期望值的左值的猜测是
But since you want the expected value between a start and +inf we need to integrate in this range, and since the PDF then no longer integrates to 1 in this smaller interval, I guess we have to normalize the result be dividing by the integral of the PDF in this range. So my guess for the left-bound expected value is
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，得到
And for the MeanResidualLife you subtract start from it, giving
MRL[start_] := expVal[start] - start
其中的情节为
Plot[MRL[start], {start, 0, 0.3}, PlotRange -> {0, All}]
看似合理，但我不是专家.所以最后我们要最小化它，即找到start，此功能是局部最小值.最小值似乎在0.05左右，但让我们从该猜测中找到一个更准确的值
Looks plausible, but I'm no expert. So finally we want to minimize it, i.e. find the start for which this function is a local minimum. The minimum seems to be around 0.05, but let's find a more exact value starting from that guess
FindMinimum[MRL[start], {start, 0.05}]
，并且在发生一些错误(您的函数未在0以下定义，所以我猜最小化器在那个禁止区域戳了一下)后，我们得到了
and after some errors (your function is not defined below 0, so I guess the minimizer pokes a little in that forbidden region) we get
{0.0418137，{开始-> 0.0584312}}
{0.0418137, {start -> 0.0584312}}
因此，最佳位置应为start = 0.0584312，平均剩余寿命为0.0418137.
So the optimum should be at start = 0.0584312 with a mean residual life of 0.0418137.
我不知道这是否正确，但这似乎是合理的.
I don't know if this is correct, but it seems plausible.
这篇关于为Mathematica中的自定义分布最小化NExpectation的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
