Ini berkaitan dengan pertanyaan sebelumnya dari bulan Juni:
Menghitung ekspektasi untuk distribusi kustom di Mathematica
Saya memiliki distribusi campuran ubahsuaian yang ditentukan menggunakan distribusi ubahsuaian kedua mengikuti garis yang dibahas @Sasha
dalam sejumlah jawaban selama setahun terakhir.
Kode yang mendefinisikan distribusi adalah sebagai berikut:
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]];
Ini memungkinkan saya untuk menyesuaikan parameter distribusi dan menghasilkan PDF dan CDF . Contoh plot:
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]
Sekarang saya telah mendefinisikan a function
untuk menghitung sisa rata-rata kehidupan (lihat pertanyaan ini untuk penjelasan).
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
Yang pertama dari ini yang tidak menetapkan batas seperti pada yang kedua membutuhkan waktu lama untuk menghitung, tetapi keduanya bekerja.
Sekarang saya perlu menemukan fungsi minimum MeanResidualLife
untuk distribusi yang sama (atau beberapa variasi) atau menguranginya.
Saya sudah mencoba beberapa variasi untuk ini:
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]
Ini tampaknya berjalan selamanya atau mengalami:
Power :: infy: Ekspresi tak terbatas 1 / 0. ditemui. >>
The MeanResidualLife
fungsi diterapkan pada sederhana tapi sama berbentuk distribusi menunjukkan bahwa ia memiliki minimal satu:
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}}]
Juga keduanya:
FindMinimum[MeanResidualLife[x, LogNormalDistribution[1.75, 0.65]], x]
FindMinimum[MeanResidualLife[x, 30, LogNormalDistribution[1.75, 0.65]], x]
beri saya jawaban (jika dengan banyak pesan terlebih dahulu) saat digunakan dengan LogNormalDistribution
.
Adakah pemikiran tentang cara membuatnya agar berfungsi untuk distribusi kustom yang dijelaskan di atas?
Apakah saya perlu menambahkan batasan atau opsi?
Apakah saya perlu mendefinisikan sesuatu yang lain dalam definisi distribusi khusus?
Mungkin FindMinimum
atau NMinimize
hanya perlu berjalan lebih lama (saya sudah menjalankannya hampir satu jam tidak berhasil). Jika demikian, apakah saya hanya perlu beberapa cara untuk mempercepat menemukan fungsi minimum? Ada saran tentang caranya?
Apakah Mathematica
ada cara lain untuk melakukan ini?
Ditambahkan 9 Feb 5:50 PM EST:
Siapa pun dapat mengunduh presentasi Oleksandr Pavlyk tentang membuat distribusi di Mathematica dari lokakarya Wolfram Technology Conference 2011 'Buat Distribusi Anda Sendiri' di sini . Unduhan termasuk notebook, 'ExampleOfParametricDistribution.nb'
yang tampaknya menjabarkan semua bagian yang diperlukan untuk membuat distribusi yang dapat digunakan seperti distribusi yang datang dengan Mathematica.
Mungkin memberikan beberapa jawaban.
Jawaban:
Sejauh yang saya lihat, masalahnya adalah (seperti yang sudah Anda tulis), yang
MeanResidualLife
membutuhkan waktu lama untuk dihitung, bahkan untuk evaluasi tunggal. Sekarang,FindMinimum
atau fungsi serupa mencoba mencari minimum untuk fungsi. Menemukan minimum memerlukan salah satu untuk menetapkan turunan pertama dari fungsi nol dan menyelesaikan solusi. Karena fungsi Anda cukup rumit (dan mungkin tidak dapat dibedakan), kemungkinan kedua adalah melakukan minimalisasi numerik, yang memerlukan banyak evaluasi fungsi Anda. Ergo, ini sangat lambat sekali.Saya sarankan untuk mencobanya tanpa sihir Mathematica.
Pertama mari kita lihat apa
MeanResidualLife
itu, seperti yang Anda definisikan.NExpectation
atauExpectation
menghitung nilai yang diharapkan . Untuk nilai yang diharapkan, kami hanya membutuhkanPDF
distribusi Anda. Mari kita ekstrak dari definisi Anda di atas ke dalam fungsi sederhana:Jika kami memplot pdf2, tampilannya persis seperti Plot Anda
Sekarang untuk nilai yang diharapkan. Jika saya memahaminya dengan benar, kita harus mengintegrasikan
x * pdf[x]
dari-inf
ke+inf
untuk nilai yang diharapkan normal.x * pdf[x]
sepertidan nilai yang diharapkan adalah
Tetapi karena Anda menginginkan nilai yang diharapkan antara a
start
dan+inf
kami perlu mengintegrasikan dalam rentang ini, dan karena PDF maka tidak lagi berintegrasi ke 1 dalam interval yang lebih kecil ini, saya kira kami harus menormalkan hasil yang akan dibagi dengan integral dari PDF di kisaran ini. Jadi tebakan saya untuk nilai yang diharapkan adalah kiriDan untuk
MeanResidualLife
Anda kurangistart
dari itu, memberiYang plot sebagai
Tampak masuk akal, tapi saya bukan ahli. Jadi akhirnya kami ingin menguranginya, yaitu menemukan
start
untuk yang fungsi ini adalah minimum lokal. Nilai minimum tampaknya sekitar 0,05, tetapi mari kita cari nilai yang lebih tepat mulai dari tebakan itudan setelah beberapa kesalahan (fungsi Anda tidak didefinisikan di bawah 0, jadi saya kira minimizer menyodok sedikit di wilayah terlarang itu) kami dapatkan
{0,0418137, {start -> 0,0584312}}
Jadi optimal harus di
start = 0.0584312
dengan sisa kehidupan rata-rata0.0418137
.Saya tidak tahu apakah ini benar, tetapi tampaknya masuk akal.
sumber