[plot] Mathematica에서 맞춤형 배포를위한 NExpectation 최소화

이것은 6 월의 초기 질문과 관련이 있습니다.

Mathematica의 사용자 정의 분포에 대한 기대 값 계산

@Sasha지난해 여러 답변에서 논의 된 내용 에 따라 두 번째 사용자 지정 배포를 사용하여 정의 된 사용자 지정 혼합 배포가 있습니다.

분포를 정의하는 코드는 다음과 같습니다.

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

이를 통해 분포 매개 변수를 맞추고 PDF 및 CDF를 생성 할 수 있습니다 . 줄거리의 예 :

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]

이제 function평균 잔류 수명을 계산 하기 위해를 정의했습니다 ( 설명 은 이 질문 을 참조하십시오 ).

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

두 번째와 같이 제한을 설정하지 않는 첫 번째는 계산하는 데 시간이 오래 걸리지 만 둘 다 작동합니다.

이제 MeanResidualLife동일한 분포 (또는 일부 변형) 에 대한 함수 의 최소값을 찾 거나 최소화해야합니다.

나는 이것에 대해 많은 변형을 시도했다.

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]

이것들은 영원히 실행되거나 다음과 같이 보입니다.

Power :: infy : 무한 표현 1 / 0이 발생했습니다. >>

MeanResidualLife함수는 하나의 최소했다는 단순한 형상이지만 마찬가지로 배신 프로그램에 적용 :

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

또한 둘 다 :

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

와 함께 사용할 때 (다양한 메시지가있는 경우 먼저) 답변을 제공하십시오 LogNormalDistribution.

위에서 설명한 사용자 정의 배포판에서이를 작동시키는 방법에 대한 생각이 있습니까?

구속 조건 또는 옵션을 추가해야합니까?

사용자 정의 배포 정의에서 다른 것을 정의해야합니까?

어쩌면 FindMinimum또는 NMinimize더 오래 실행해야 할 수도 있습니다 (아무런 시간도 거의 없어야합니다). 그렇다면 함수의 최소값을 찾는 속도를 높이는 방법이 필요합니까? 방법에 대한 제안?

않습니다 Mathematica이 작업을 수행하는 또 다른 방법이 있나요?

추가됨 2 월 9 일 5:50 PM EST :

누구나 Wolfram Technology Conference 2011 워크숍 ‘Create Your Own Distribution’에서 Mathematica에서 배포판 작성에 관한 Oleksandr Pavlyk의 프리젠 테이션을 다운로드 할 수 있습니다 . 다운로드에는 노트북이 포함되어 있는데, 'ExampleOfParametricDistribution.nb'Mathematica와 함께 제공되는 배포판처럼 사용할 수있는 배포판을 만드는 데 필요한 모든 부분을 레이아웃하는 것 같습니다.

일부 답변이 제공 될 수 있습니다.

답변

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;

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

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

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

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

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

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

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