Minimización NExpectation para una distribución personalizada en Mathematica

Question

Esto se refiere a una pregunta anterior de la parte posterior en junio:

Calculating expectation for a custom distribution in Mathematica

I tiene una distribución de mezclado de encargo define utilizando una segunda distribución personalizada siguiente a lo largo de las líneas discutidas por @Sasha en varias respuestas durante el año pasado.

Código definir las distribuciones de la siguiente manera:

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

Esto me permite adaptarse a los parámetros de distribución y generar de PDF y CDF de. Un ejemplo de las parcelas:

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] 

Ahora he definido un function para calcular la vida media residual (ver this question para una explicación).

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 

El primero de ellos que no establece un límite al igual que en la segunda toma mucho tiempo para calcular, pero que ambos trabajen.

Ahora necesito encontrar el mínimo de la función MeanResidualLife para la misma distribución (o alguna variación de la misma) o minimizarla.

He intentado una serie de variaciones sobre esto:

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] 

Estos tampoco parecen funcionar para siempre o correr en:

Poder :: INFY: expresión Infinito 1/0. encontrado . >>

La función MeanResidualLife aplica a un simple pero la distribución de forma similar muestra que tiene un único mínimo:

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

también tanto:

FindMinimum[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], x] 
FindMinimum[MeanResidualLife[x, 30, LogNormalDistribution[1.75, 0.65]], x] 

dar respuestas (si con un montón de mensajes primero) cuando se usa con el LogNormalDistribution.

¿Alguna idea sobre cómo hacer que esto funcione para la distribución personalizada descrita anteriormente?

¿Debo agregar restricciones u opciones?

¿Debo definir algo más en las definiciones de las distribuciones personalizadas?

Tal vez el FindMinimum o el NMinimize solo necesiten ejecutarse más (los he ejecutado durante casi una hora en vano). Si es así, ¿necesito alguna forma de acelerar la búsqueda del mínimo de la función? ¿Alguna sugerencia sobre cómo?

¿Tiene Mathematica otra forma de hacerlo?

añadieron 9 Feb 17:50 EST:

Cualquier persona puede descargar la presentación de Oleksandr Pavlyk sobre la creación de distribuciones en Mathematica de la Conferencia de Tecnología Wolfram 2011 Taller 'Su distribución' here. Las descargas incluyen el cuaderno, 'ExampleOfParametricDistribution.nb' que parece presentar todas las piezas necesarias para crear una distribución que se pueda usar como las distribuciones que vienen con Mathematica.

Puede proporcionar parte de la respuesta.

Answer 1

Por lo que veo, el problema es (como usted ya escribió), que MeanResidualLife tarda mucho tiempo en computarse, incluso para una evaluación única. Ahora, el FindMinimum o funciones similares intentan encontrar un mínimo a la función. Encontrar un mínimo requiere establecer la primera derivada de la función cero y resolver una solución. Como su función es bastante complicada (y probablemente no diferenciable), la segunda posibilidad es realizar una minimización numérica, que requiere muchas evaluaciones de su función. Ergo, es muy, muy lento.

Sugeriría probar sin la magia de Mathematica.

Primero veamos qué es MeanResidualLife, como usted lo definió. NExpectation o Expectation calcule el expected value. Para el valor esperado, solo necesitamos el PDF de su distribución. Vamos a extraerlo de su definición anterior en funciones simples:

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; 

Si dibujamos pdf2 se ve exactamente como su casa,

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

Ahora al valor esperado. Si lo entiendo correctamente, tenemos que integrar x * pdf[x] de -inf a +inf para un valor esperado normal.

x * pdf[x] parece

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

y el valor esperado es

NIntegrate[pdf2[3.77, 1.34, -2.65, 0.40, x]*x, {x, 0, \[Infinity]}] 
Out= 0.0596504 

Pero ya que desea que el valor esperado entre un start y +inf tenemos que integrar en este rango , y dado que el PDF ya no se integra a 1 en este intervalo más pequeño, supongo que tenemos e para normalizar el resultado dividiéndolo por la integral del PDF en este rango.Así que yo creo para el valor esperado izquierda unida es

Y para el MeanResidualLife se resta start de ella, dando

MRL[start_] := expVal[start] - start 

que traza como

Plot[MRL[start], {start, 0, 0.3}, PlotRange -> {0, All}] 

Parece plausible, pero no soy un experto. Entonces, finalmente queremos minimizarlo, es decir, encontrar el start para el cual esta función es un mínimo local. El mínimo parece ser alrededor de 0,05, pero vamos a encontrar un valor más exacto a partir de esa conjetura

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

y después de algunos errores (su función no está definida por debajo de 0, así que supongo que el minimizador asoma un poco en ese prohibida región) obtenemos

{0,0418137, {Inicio -> 0.0584312}}

lo tanto el óptimo debe ser de al start = 0.0584312 con una vida media residual de 0.0418137.

No sé si esto es correcto, pero parece plausible.