plot wolfram-mathematica probability calculus mathematica-8

plot - Minimizando la expectativa de una distribución personalizada en Mathematica



wolfram-mathematica probability (1)

Por lo que veo, el problema es (como ya escribiste), que MeanResidualLife tarda mucho tiempo en computarse, incluso para una sola evaluación. Ahora, el FindMinimum o funciones 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. Debido a que su función es bastante complicada (y probablemente no sea diferenciable), la segunda posibilidad es realizar una minimización numérica, lo que requiere muchas evaluaciones de su función. Ergo, es muy muy lento.

Te sugiero que lo pruebes sin la magia de Mathematica.

Primero veamos qué es el MeanResidualLife , como lo MeanResidualLife . NExpectation Expectation o la Expectation calcula el valor esperado . Para el valor esperado, solo necesitamos el PDF de su distribución. Extraigámoslo 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 parcela

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

Ahora al valor esperado. Si lo comprendo correctamente, debemos integrar x * pdf[x] de -inf a +inf para obtener un valor normal esperado.

x * pdf[x] parece a

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 como desea que el valor esperado esté entre un start y un +inf , necesitamos integrarlo en este rango, y dado que el PDF ya no se integra a 1 en este intervalo más pequeño, creo que tenemos que normalizar el resultado dividiéndolo por la integral de El PDF en este rango. Así que mi conjetura para el valor esperado límite a la izquierda 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 MeanResidualLife le resta start de él, dando

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

Que parcelas como

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

Parece plausible, pero no soy un experto. Finalmente, queremos minimizarlo, es decir, encontrar el start para el cual 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 por debajo de 0, así que supongo que el minimizador se asoma un poco en esa región prohibida) obtenemos

{0.0418137, {inicio -> 0.0584312}}

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

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

Esto se relaciona con una pregunta anterior de junio:

Cálculo de la expectativa de una distribución personalizada en Mathematica.

Tengo una distribución mixta personalizada definida utilizando una segunda distribución personalizada que sigue las líneas comentadas por @Sasha en 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 archivos PDF y CDF . 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 una function para calcular la vida residual media (consulte 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 función MeanResidualLife para la misma distribución (o alguna variación de ella) o minimizarla.

He intentado una serie de variaciones en 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 o bien parecen funcionar para siempre o se encuentran con:

Potencia :: infy: Expresión infinita 1 / 0. encontrado. >>

La función MeanResidualLife aplicada a una distribución más simple pero con una 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}}]

También ambos:

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

dame respuestas (si con un montón de mensajes primero) cuando se usa con 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 FindMinimum o NMinimize solo necesite correr más tiempo (los he ejecutado durante casi una hora sin éxito). Si es así, ¿necesito alguna forma de acelerar la búsqueda del mínimo de la función? ¿Alguna sugerencia sobre cómo?

¿ Mathematica tiene otra forma de hacer esto?

Añadido el 9 de febrero 5:50 p.m. EST:

Cualquiera puede descargar la presentación de Oleksandr Pavlyk sobre la creación de distribuciones en Mathematica desde el taller de la Conferencia de Tecnología Wolfram 2011 ''Crea tu propia distribución'' here . Las descargas incluyen el cuaderno, ''ExampleOfParametricDistribution.nb'' que parece ''ExampleOfParametricDistribution.nb'' todas las piezas necesarias para crear una distribución que se puede usar como las distribuciones que vienen con Mathematica.

Puede proporcionar algo de la respuesta.