Dies bezieht sich auf eine frühere Frage aus dem Juni:
Berechnung der Erwartung für eine benutzerdefinierte Verteilung in Mathematica
Ich habe eine benutzerdefinierte gemischte Verteilung definiert, die unter Verwendung einer zweiten benutzerdefinierten Verteilung definiert wird, die den @Sasha
in einer Reihe von Antworten des letzten Jahres erörterten Linien folgt .
Der Code, der die Verteilungen definiert, lautet wie folgt:
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]];
Dadurch kann ich Verteilungsparameter anpassen und PDFs und CDFs generieren . Ein Beispiel für die Handlungen:
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]
Jetzt habe ich a definiert function
, um die mittlere Restlebensdauer zu berechnen ( eine Erklärung finden Sie in dieser Frage ).
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
Die erste von diesen, die kein Limit wie in der zweiten festlegt, dauert lange, aber beide funktionieren.
Jetzt muss ich das Minimum der MeanResidualLife
Funktion für dieselbe Verteilung (oder eine Variation davon) finden oder minimieren.
Ich habe eine Reihe von Variationen davon ausprobiert:
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]
Diese scheinen entweder für immer zu laufen oder treffen auf:
Power :: infy: Unendlicher Ausdruck 1 / 0. angetroffen. >>
Die MeanResidualLife
Funktion, die auf eine einfachere, aber ähnlich geformte Verteilung angewendet wird, zeigt, dass sie ein einziges Minimum hat:
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}}]
Auch beides:
FindMinimum[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], x]
FindMinimum[MeanResidualLife[x, 30, LogNormalDistribution[1.75, 0.65]], x]
Geben Sie mir Antworten (wenn zuerst mit einer Reihe von Nachrichten), wenn Sie mit dem verwendet werden LogNormalDistribution
.
Irgendwelche Gedanken darüber, wie dies für die oben beschriebene benutzerdefinierte Distribution funktioniert?
Muss ich Einschränkungen oder Optionen hinzufügen?
Muss ich in den Definitionen der benutzerdefinierten Distributionen etwas anderes definieren?
Vielleicht müssen die FindMinimum
oder NMinimize
müssen nur länger laufen (ich habe sie fast eine Stunde ohne Erfolg laufen lassen). Wenn ja, brauche ich nur eine Möglichkeit, um das Minimum der Funktion zu finden? Irgendwelche Vorschläge wie?
Hat Mathematica
eine andere Möglichkeit, dies zu tun?
Hinzugefügt am 9. Februar, 17:50 Uhr EST:
Jeder kann die Präsentation von Oleksandr Pavlyk über das Erstellen von Distributionen in Mathematica vom Workshop 'Create Your Own Distribution' der Wolfram Technology Conference 2011 hier herunterladen . Zu den Downloads gehört das Notizbuch, 'ExampleOfParametricDistribution.nb'
das alle Teile enthält, die zum Erstellen einer Distribution erforderlich sind, die wie die mit Mathematica gelieferten Distributionen verwendet werden kann.
Es kann einen Teil der Antwort liefern.
Antworten:
Soweit ich sehe, besteht das Problem (wie Sie bereits geschrieben haben) darin, dass
MeanResidualLife
die Berechnung selbst für eine einzelne Auswertung lange dauert. NunFindMinimum
versuchen die oder ähnliche Funktionen, ein Minimum für die Funktion zu finden. Um ein Minimum zu finden, muss entweder die erste Ableitung der Funktion Null gesetzt und nach einer Lösung gesucht werden. Da Ihre Funktion ziemlich kompliziert (und wahrscheinlich nicht differenzierbar) ist, besteht die zweite Möglichkeit darin, eine numerische Minimierung durchzuführen, die viele Bewertungen Ihrer Funktion erfordert. Ergo ist es sehr sehr langsam.Ich würde vorschlagen, es ohne Mathematica-Magie zu versuchen.
Lassen Sie uns zuerst sehen, was das
MeanResidualLife
ist, wie Sie es definiert haben.NExpectation
oderExpectation
berechnen Sie den erwarteten Wert . Für den erwarteten Wert benötigen wir nur diePDF
Ihrer Distribution. Lassen Sie es uns aus Ihrer obigen Definition in einfache Funktionen extrahieren:Wenn wir pdf2 zeichnen, sieht es genauso aus wie Ihr Plot
Nun zum erwarteten Wert. Wenn ich es richtig verstehe, müssen wir für einen normalen Erwartungswert
x * pdf[x]
von-inf
bis integrieren+inf
.x * pdf[x]
sieht aus wieund der erwartete Wert ist
Da Sie jedoch den erwarteten Wert zwischen a
start
und+inf
in diesen Bereich integrieren möchten und das PDF in diesem kleineren Intervall nicht mehr auf 1 integriert werden muss, müssen wir das Ergebnis normalisieren, indem wir es durch das Integral des PDF in dividieren dieser Bereich. Meine Vermutung für den linksgebundenen erwarteten Wert ist alsoUnd für die
MeanResidualLife
du davon subtrahierststart
, gibstWelche Handlungen als
Sieht plausibel aus, aber ich bin kein Experte. Schließlich wollen wir es minimieren, dh das finden,
start
für das diese Funktion ein lokales Minimum ist. Das Minimum scheint bei 0,05 zu liegen, aber lassen Sie uns ausgehend von dieser Vermutung einen genaueren Wert findenund nach einigen Fehlern (Ihre Funktion ist nicht unter 0 definiert, also denke ich, dass der Minimierer ein wenig in dieser verbotenen Region steckt) bekommen wir
{0.0418137, {start -> 0.0584312}}
Das Optimum sollte also bei
start = 0.0584312
einer mittleren Restlebensdauer von liegen0.0418137
.Ich weiß nicht, ob das richtig ist, aber es scheint plausibel.
quelle