Minimizando NExpectation para una distribución personalizada en Mathematica

238

Esto se relaciona con una pregunta anterior de junio:

Cálculo de expectativas para una distribución personalizada en Mathematica

Tengo una distribución mixta personalizada definida usando una segunda distribución personalizada siguiendo las líneas discutidas @Sashaen varias respuestas durante el año pasado.

El código que define las distribuciones sigue:

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 ajustar los parámetros de distribución y generar PDF y CDF . Un ejemplo de las tramas:

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]

ingrese la descripción de la imagen aquí

Ahora he definido a functionpara calcular la vida residual media (vea esta pregunta para obtener 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 estos que no establece un límite como en el segundo toma mucho tiempo para calcular, pero ambos funcionan.

Ahora necesito encontrar el mínimo de la MeanResidualLifefunción 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 parecen ejecutarse para siempre o encontrarse con:

Poder :: infy: Expresión infinita 1 / 0. encontrada. >>

La MeanResidualLifefunción aplicada a una distribución más simple pero de forma similar muestra que tiene un mínimo único:

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

ingrese la descripción de la imagen aquí

También ambos:

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

darme respuestas (si es 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?

¿Necesito agregar restricciones u opciones?

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

Tal vez el FindMinimumo NMinimizesimplemente necesito correr más tiempo (los he ejecutado 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?

¿ MathematicaTiene otra forma de hacer esto?

Añadido 9 feb 5:50 PM EST:

Cualquiera puede descargar la presentación de Oleksandr Pavlyk sobre la creación de distribuciones en Mathematica desde el taller Wolfram Technology Conference 2011 'Cree su propia distribución' aquí . Las descargas incluyen el cuaderno, 'ExampleOfParametricDistribution.nb'que parece presentar todas las piezas necesarias para crear una distribución que se puede usar como las distribuciones que vienen con Mathematica.

Puede proporcionar algunas de las respuestas.

plot wolfram-mathematica probability calculus mathematica-8 Jagra
fuente
99
No es experto en Mathematica, pero he encontrado problemas similares en otros lugares. Parece que tiene problemas cuando su dominio comienza en 0. Intente comenzar en 0.1 y más y vea qué sucede.
Makketronix
77
@ Makketronix - Gracias por esto. Divertida sincronía, dado que comencé a revisar esto después de 3 años.
Jagra
8
No estoy seguro de poder ayudarlo, pero podría intentar preguntar en el stackoverflow específico de Mathematica . ¡La mejor de las suertes!
Olivia Stork
44
¿Intentó: reference.wolfram.com/language/ref/Expectation.html ?
Cplusplusplus el
1
Hay muchos artículos al respecto en zbmath.org Búsqueda de expectativas
Ivan V

Respuestas:

11

Por lo que veo, el problema es (como ya escribió), que MeanResidualLifelleva mucho tiempo calcular, incluso para una sola evaluación. Ahora, las FindMinimumfunciones o similares intentan encontrar un mínimo para la función. Encontrar un mínimo requiere establecer la primera derivada de la función cero y resolver una solución. Dado que su función es bastante complicada (y probablemente no sea diferenciable), la segunda posibilidad es hacer una minimización numérica, que requiere muchas evaluaciones de su función. Ergo, es muy muy lento.

Sugeriría probarlo sin la magia de Mathematica.

Primero veamos qué MeanResidualLifees, como lo definiste. NExpectationo Expectationcalcular el valor esperado . Para el valor esperado, solo necesitamos el PDFde 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 trazamos pdf2 se ve exactamente como su trama

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

Parcela de PDF

Ahora al valor esperado. Si lo entiendo correctamente, tenemos que integrar x * pdf[x]de -infa +infpara obtener un valor normal esperado.

x * pdf[x] parece

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

Parcela de x * PDF

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 una starty +infnecesitamos integrar en este rango, y puesto que el PDF se integra a continuación, ya no a 1 en este intervalo más pequeño, supongo que tenemos para normalizar el resultado se divide por la integral de la PDF en este rango. Entonces, mi suposición para el valor esperado del límite izquierdo es

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

Y para el MeanResidualLifeque restas startde él, dando

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

Que trama como

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

Trama de vida residual media

Parece plausible, pero no soy un experto. Así que finalmente queremos minimizarlo, es decir, encontrar startpara qué esta función es un mínimo local. El mínimo parece estar alrededor de 0.05, pero encontremos un valor más exacto a partir de esa suposición

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

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

{0.0418137, {inicio -> 0.0584312}}

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

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

azt
fuente
+1: acabo de ver esto, así que tendré que resolverlo, pero creo que la forma en que dividiste el problema en pasos solucionables tiene mucho sentido. Además, la trama de su función MRL, sin duda, parece acertada. Muchas gracias, volveré sobre esto tan pronto como pueda hacer tiempo para estudiar su respuesta.
Jagra