在 Mathematica 中最小化自定义分布的 NExpectation

Posted 2023-02-16

【中文标题】在 Mathematica 中最小化自定义分布的 NExpectation

【英文标题】：Minimizing NExpectation for a custom distribution in Mathematica

【发布时间】：2012-03-01 10:21:56

【问题描述】：

这与 6 月份的一个较早的问题有关：

Calculating expectation for a custom distribution in 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 来计算平均剩余寿命（请参阅this question 以获得解释）。

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：

任何人都可以从 Wolfram 技术会议 2011 研讨会“创建您自己的发行版”here 下载 Oleksandr Pavlyk 的关于在 Mathematica 中创建发行版的演示文稿。下载内容包括笔记本 'ExampleOfParametricDistribution.nb'，它似乎列出了创建分发所需的所有部分，可以像 Mathematica 附带的分发一样使用。

它可能会提供一些答案。

【问题讨论】：

不是 Mathematica 专家，但我在其他地方也遇到过类似的问题。当您的域从 0 开始时，您似乎遇到了问题。尝试从 0.1 及更高版本开始，看看会发生什么。

@Makketronix -- 谢谢你。有趣的同步性，因为我在 3 年后开始重新审视它。

我不确定能否帮助您，但您可以尝试通过Mathematica-specific *** 询问。祝你好运！

你试过了吗：reference.wolfram.com/language/ref/Expectation.html？

zbmath.org上有一堆关于它的文章搜索期待

【参考方案1】：

据我所知，问题是（正如您已经写过的），MeanResidualLife 需要很长时间来计算，即使是单次评估也是如此。现在，FindMinimum 或类似函数试图找到该函数的最小值。找到最小值需要将函数的一阶导数设置为零并求解。由于您的函数非常复杂（并且可能不可微），第二种可能性是进行数值最小化，这需要对您的函数进行多次评估。因此，它非常非常慢。

我建议在没有 Mathematica 魔法的情况下尝试一下。

首先让我们看看你定义的MeanResidualLife 是什么。 NExpectation 或 Expectation 计算 expected value。 对于期望值，我们只需要您分配的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, 开始 -> 0.0584312

因此最佳值应为start = 0.0584312，平均剩余寿命为0.0418137。

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

【讨论】：

+1 -- 刚看到这个，所以我需要解决它，但我认为你将问题划分为可解决步骤的方式很有意义。此外，您的 MRL 函数图看起来很准确。非常感谢，我会尽快抽出时间研究您的答案。