Minimizzare NExpectation per una distribuzione personalizzato in Mathematica

Question

sopra si riferisce a una domanda precedente da a giugno:

Calculating expectation for a custom distribution in Mathematica

Ho una distribuzione mista personalizzato definito utilizzando una seconda distribuzione personalizzata seguendo lungo le linee discussi da @Sasha in un certo numero di risposte nell'ultimo anno.

codice che definisce le distribuzioni segue:

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

Questo mi permette di adattare i parametri di distribuzione e generare di PDF e CDF di. Un esempio delle piazzole:

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] 

Ora ho definito un function per calcolare la vita residua media (vedi this question per una spiegazione).

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 

Il primo di questi che non fissa un limite nella seconda richiede molto tempo per calcolare, ma entrambi lavorano.

Ora ho bisogno di trovare il minimo della funzione MeanResidualLife per la stessa distribuzione (o qualche variazione di esso) o ridurla a icona.

Ho provato un certo numero di variazioni su questo:

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] 

Questi due sembrano correre per sempre o incorrere in:

Potenza :: INFY: espressione Infinite 1/0. incontrato . >>

La funzione MeanResidualLife applicata ad un semplice ma la distribuzione di forma simile dimostra che ha un singolo minimo:

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

anche entrambi:

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

give Io rispondo (se con un gruppo di messaggi prima) quando usato con lo LogNormalDistribution.

Qualche idea su come farlo funzionare per la distribuzione personalizzata descritta sopra?

Devo aggiungere vincoli o opzioni?

Devo definire qualcos'altro nelle definizioni delle distribuzioni personalizzate?

Forse il FindMinimum o NMinimize devono solo essere più lunghi (li ho eseguiti quasi un'ora senza risultati). In tal caso, ho solo bisogno di un modo per accelerare la ricerca del minimo della funzione? Qualche suggerimento su come?

Mathematica ha un altro modo per farlo?

Aggiunto 9 Feb 17:50 EST:

chiunque può scaricare la presentazione di Oleksandr Pavlyk di creare distribuzioni in Mathematica dalla Wolfram Technology Conference 2011 workshop 'creare la vostra distribuzione' here. I download includono il notebook, 'ExampleOfParametricDistribution.nb', che sembra disporre di tutti i pezzi necessari per creare una distribuzione utilizzabile come le distribuzioni fornite con Mathematica.

Può fornire una risposta.

Answer 1

Per quanto vedo, il problema è (come hai già scritto), che MeanResidualLife richiede molto tempo per il calcolo, anche per una singola valutazione. Ora, le funzioni FindMinimum o simili cercano di trovare un minimo per la funzione. Per trovare un minimo è necessario impostare la prima derivata della funzione zero e risolvere una soluzione. Dato che la tua funzione è piuttosto complicata (e probabilmente non è differenziabile), la seconda possibilità è quella di eseguire una minimizzazione numerica, che richiede molte valutazioni della tua funzione. Ergo, è molto molto lento.

Suggerirei di provarlo senza la magia di Mathematica.

Per prima cosa vediamo cosa è lo MeanResidualLife, come lo hai definito. NExpectation o Expectation calcolare lo expected value. Per il valore atteso, abbiamo solo bisogno dello PDF della vostra distribuzione. Diamo estrarla dalla vostra definizione di cui sopra in funzioni semplici:

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; 

Se tracciamo pdf2 sembra esattamente le Trama

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

Ora al valore atteso. Se ho capito bene, dobbiamo integrare x * pdf[x] da -inf a +inf per un normale valore atteso.

x * pdf[x] sembra

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

e il valore atteso è

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

Ma dal momento che si desidera che il valore atteso tra un start e +inf abbiamo bisogno di integrare in questa fascia e dal momento che il PDF non si integra più a 1 in questo intervallo più piccolo, suppongo che abbiamo e per normalizzare il risultato dividere per l'integrale del PDF in questo intervallo.Quindi la mia ipotesi per il valore atteso rilegato a sinistra è

E per la MeanResidualLife si sottrae start da esso, dando

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

che riporta come

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

Sembra plausibile, ma non sono esperto. Quindi, alla fine, vogliamo ridurlo al minimo, vale a dire trovare lo start per il quale questa funzione è un minimo locale. Il minimo sembra essere intorno a 0,05, ma cerchiamo di trovare un valore più preciso a partire da tale supposizione

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

e dopo alcuni errori (la funzione non è definita sotto 0, quindi credo che il Minimizer sporge un po 'in quella proibita regione) otteniamo

{,0418,137 mila, {Avvio -> 0,0584,312 mila}}

Così l'ottimale dovrebbe essere a start = 0.0584312 con una vita residua media di 0.0418137.

Non so se questo è corretto, ma sembra plausibile.