Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: More Inquiries

  • To: mathgroup at smc.vnet.net
  • Subject: [mg91250] Re: [mg91163] More Inquiries
  • From: Bob Hanlon <hanlonr at cox.net>
  • Date: Mon, 11 Aug 2008 06:07:45 -0400 (EDT)
  • Reply-to: hanlonr at cox.net

Clear[BurrDistribution];

BurrDistribution::usage = 
  "BurrDistribution[c, k] is continuous probability distribution for a \
non-negative random variable.";

BurrDistribution /: PDF[
   BurrDistribution[c_, k_], x_] :=
  (c*k)*(x^(c - 1)/(1 + x^c)^(k + 1));

Verify density:

Integrate[PDF[BurrDistribution[c, k], x], {x, 0, Infinity},
 Assumptions -> {c > 0, k > 0}]

1

Mean

Integrate[x*PDF[BurrDistribution[c, k], x], {x, 0, Infinity},
 Assumptions -> {c > 0, k > 0, c*k > 1}]

(Gamma[1 + 1/c]*Gamma[k - 1/c])/Gamma[k]

BurrDistribution /: Mean[
   BurrDistribution[c_, k_]] :=
  k*Beta[1 + 1/c, k - 1/c];

Mode

Simplify[Reduce[{D[PDF[BurrDistribution[c, k], x], x] == 0, c > 0, k > 0}, x, 
  Reals], {k > 0, x > 0}]

c > 1 && ((c - 1)/(c*k + 1))^(1/c) == x

BurrDistribution /: Mode[
   BurrDistribution[c_, k_]] :=
  ((c - 1)/(c*k + 1))^(1/c);

CDF

Integrate[PDF[BurrDistribution[c, k], t], {t, 0, x},
 Assumptions -> {c > 0, k > 0}]

1 - (x^c + 1)^(-k)

BurrDistribution /: CDF[
   BurrDistribution[c_, k_], x_] :=
  1 - (x^c + 1)^(-k);

Limit[CDF[BurrDistribution[c, k], x], x -> 0,
 Assumptions -> {c > 0, k > 0}]

0

Limit[CDF[BurrDistribution[c, k], x], x -> Infinity,
 Assumptions -> {c > 0, k > 0}]

1

Inverse CDF

Off[Solve::ifun];

x /. Solve[CDF[BurrDistribution[c, k], x] == q, x][[1]]

((1 - q)^(-k^(-1)) - 1)^(1/c)

BurrDistribution /: InverseCDF[
      BurrDistribution[c_, k_], q_] :=
    ((1 - q)^(-1/k) - 1)^(1/c);

Median

InverseCDF[BurrDistribution[c, k], 1/2]

(-1 + 2^(1/k))^(1/c)

BurrDistribution /: Median[
   BurrDistribution[c_, k_]] :=
  (2^(1/k) - 1)^(1/c);

RandomReal

BurrDistribution /: RandomReal[
   BurrDistribution[c_?NumericQ, k_?NumericQ]] :=
  
  InverseCDF[BurrDistribution[c, k], RandomReal[]];

BurrDistribution /: RandomReal[
   BurrDistribution[c_?NumericQ, k_?NumericQ], 
   n_Integer?Positive] :=
  
  Table[InverseCDF[BurrDistribution[c, k], RandomReal[]], {n}];

Parameter Estimation

x2 = {5, 6, 5.5, 5, 8.4, 4.5, 4, 6, 7, 6, 7, 8, 5.4, 5, 4, 3.5, 4.5,
   6.5, 5, 3.5, 5.7, 5.8, 5, 4.8, 4.5, 5.8, 3.2, 6, 4, 6, 7, 6.5,
   10.7, 7, 4.7, 7, 8.3, 10, 6.5, 4.9, 3.4, 9, 6, 3.1, 5, 4.8, 4,
   6, 6, 5.6, 4.2, 4.3, 4.5, 8.4, 8.6, 5.8, 6.8, 4.8, 3.4, 8, 8,
   8.3, 8, 4.5, 6, 4.5, 7, 7, 8};

Your data appears to have an offset from the origin. Perhaps you should model an offset.

{xmin = Min[x2], xmean = Mean[x2], xmax = Max[x2]}

{3.1,5.88696,10.7}

llf = Total[Log[PDF[BurrDistribution[c, k], #] & /@ (x2 - d)]];

paramEst = NMinimize[{-llf,
    c > 0 && k > 0 && c*k > 1 && 
     0 <= d < xmin &&
     (Mean[BurrDistribution[c, k]] + d == xmean)}, {c, 
    k, d}][[2]]

{c->2.59211,d->3.02113,k->0.568387}

I have used the mean to reduce one of the degrees of freedom and stabilize the fit.

Needs["Histograms`"];

Show[{Histogram[x2, HistogramScale -> 1],
  Plot[PDF[BurrDistribution[c, k], (x - d)] /. paramEst,
   {x, 0, Ceiling[xmax]},
   PlotStyle -> Red, PlotRange -> All]}]

The Burr distribution does not appear to fit the data very well. Since the Burr distribution clusters near the origin perhaps it would be a better match to the log of the data.

logx2 = Log[x2];

{xmin = Min[logx2], xmean = Mean[logx2], xmax = Max[logx2]}

{1.1314,1.73226,2.37024}

llf = Total[Log[PDF[BurrDistribution[c, k], #] & /@ (logx2 - d)]];

paramEst = NMinimize[{-llf,
    c > 0 && k > 0 && c*k > 1 && 
     0 <= d < xmin &&
     (Mean[BurrDistribution[c, k]] + d == xmean)}, {c, 
    k, d}][[2]]

{c->3.70856,d->0.937691,k->2.21586}

Show[{Histogram[logx2, HistogramScale -> 1],
  Plot[PDF[BurrDistribution[c, k], (x - d)] /. paramEst,
   {x, 0, Ceiling[xmax]},
   PlotStyle -> Red, PlotRange -> All]}]

The Burr distribution is a better fit to the log of your data.

x1 = {3, 5, 6, 14, 11, 7, 7, 10, 11, 5, 4, 19, 9, 2, 8, 5, 6, 6, 5,
   5, 4, 5, 5, 6, 8, 5, 8, 6, 5, 16, 5, 18, 16, 21, 6, 3, 16, 8,
   3, 8, 11, 2, 3, 4, 8, 7, 9, 10, 8, 11, 8, 10, 9, 12, 12, 9, 6,
   12, 3, 9, 14, 7, 4, 13, 8, 14, 5, 8, 2};

{xmin = Min[x1], xmean = Mean[x1], xmax = Max[x1]} // N

{2.,8.08696,21.}

llf = Total[Log[PDF[BurrDistribution[c, k], #] & /@ (x1 - d)]];

paramEst = NMinimize[{-llf,
    c > 0 && k > 0 && c*k > 1 && 
     0 <= d < xmin &&
     (Mean[BurrDistribution[c, k]] + d == xmean)}, {c, 
    k, d}][[2]]

{c->2.13746,d->1.86989,k->0.552643}

Show[{Histogram[x1, HistogramScale -> 1],
  Plot[PDF[BurrDistribution[c, k], (x - d)] /. paramEst,
   {x, 0, Ceiling[xmax]},
   PlotStyle -> Red, PlotRange -> All]}]

Again, the Burr distribution does not appear to fit the data very well. Look at the log of the data.

logx1 = Log[{3, 5, 6, 14, 11, 7, 7, 10, 11, 5, 4, 19, 9, 2, 8, 5, 6, 6, 5,
    5, 4, 5, 5, 6, 8, 5, 8, 6, 5, 16, 5, 18, 16, 21, 6, 3, 16, 8,
    3, 8, 11, 2, 3, 4, 8, 7, 9, 10, 8, 11, 8, 10, 9, 12, 12, 9, 6,
    12, 3, 9, 14, 7, 4, 13, 8, 14, 5, 8, 2}];

{xmin = Min[logx1], xmean = Mean[logx1], xmax = Max[logx1]} // N

{0.693147,1.95016,3.04452}

llf = Total[Log[PDF[BurrDistribution[c, k], #] & /@ (logx1 - d)]];

paramEst = NMinimize[{-llf,
    c > 0 && k > 0 && c*k > 1 && 
     0 <= d < xmin &&
     (Mean[BurrDistribution[c, k]] + d == xmean)}, {c, 
    k, d}][[2]]

{c->3.55713,d->0.551778,k->0.742938}

Show[{Histogram[logx1, HistogramScale -> 1],
  Plot[PDF[BurrDistribution[c, k], (x - d)] /. paramEst,
   {x, 0, Ceiling[xmax]},
   PlotStyle -> Red, PlotRange -> All]}]


Bob Hanlon

---- Dewi Anggraini <dewi_anggraini at student.rmit.edu.au> wrote: 

=============
Hi All

I've tried this following formula to Gamma distribution, it estimates the unknown parameters very well

n = 69;
x1 = {3, 5, 6, 14, 11, 7, 7, 10, 11, 5, 4, 19, 9, 2, 8, 5, 6, 6, 5,
   5,
       4, 5, 5, 6, 8, 5, 8, 6, 5, 16, 5, 18, 16, 21, 6, 3, 16, 8, 3,
   8,
       11, 2, 3, 4, 8, 7, 9, 10, 8, 11, 8, 10, 9, 12, 12, 9, 6, 12,
   3, 9,
       14, 7, 4, 13, 8, 14, 5, 8, 2};
pdf = PDF[GammaDistribution[\[Lambda], \[Beta]], x1]
logl = Plus @@ Log[pdf]
maxlogl = FindMinimum[-logl, {\[Lambda], 1}, {\[Beta], 1}]
mle = maxlogl[[2]]

{190.145, {\[Lambda] -> 3.72763, \[Beta] -> 2.16947}}

{\[Lambda] -> 3.72763, \[Beta] -> 2.16947}


and also the program produces the algebraic form of Gamma distribution by the following command and it calculates the probability of non-conformance (above the specification limit). The program also can produce the pdf graph.


PDF[GammaDistribution[\[Lambda], \[Beta]], t]

(\[ExponentialE]^-t/\[Beta] t^(-1 + \[Lambda]) \
\[Beta]^-\[Lambda])/Gamma[\[Lambda]]

Integrate[
 PDF[GammaDistribution[3.7276258602234646, 2.169465718015002], t], {t,
   8, Infinity}]

0.439226

Plot[PDF[GammaDistribution[3.7276258602234646, 2.169465718015002],
  t], {t, 2, 21}]


However, when I apply the same command towards Burr distribution (the following program below), especially for the case of producing algebraic formula of the pdf, calculating the probability of non-conformance and drawing the pdf graph (the last three commands), the program does not work very well. Do I have done something wrong in the program?



n = 69;
x1 = {3, 5, 6, 14, 11, 7, 7, 10, 11, 5, 4, 19, 9, 2, 8, 5, 6, 6, 5,
   5,
       4, 5, 5, 6, 8, 5, 8, 6, 5, 16, 5, 18, 16, 21, 6, 3, 16, 8, 3,
   8,
       11, 2, 3, 4, 8, 7, 9, 10, 8, 11, 8, 10, 9, 12, 12, 9, 6, 12,
   3, 9,
       14, 7, 4, 13, 8, 14, 5, 8, 2};
BurrDistribution[x1_, c_,
  k_] := (c*k)*(x1^(c - 1)/(1 + x1^c)^(k + 1))
pdf = BurrDistribution[x1, c, k]
logl = Plus @@ Log[pdf]
maxlogl = FindMinimum[{-logl, c > 0 && k > 0}, {c, 1}, {k, 2}, 
     MaxIterations -> 1000]
mle = maxlogl[[2]]

{249.647, {c -> 37.8115, k -> 0.0135614}}

{c -> 37.8115, k -> 0.0135614}


PDF[BurrDistribution[c, k], t]

Integrate[
 PDF[BurrDistribution[37.81151579009424, 0.01356141249769735], t] {t, 
   21, Infinity}]

Plot[PDF[BurrDistribution[37.81151579009424, 0.01356141249769735], 
  t], {t, 2, 21}]


I already got assistance from some of you about gaining MLE of Burr distribution for my data and as recommended now I can run it very well.

Now, please assist me again to find where I got wrong in my second case here (to find the algebraic form of Burr, do the integration and draw the graph). Thank You.


Kindly Regards,
Dewi






  • Prev by Date: RE: Re: sum question
  • Next by Date: Re: Solving nonlinear inequality constraints
  • Previous by thread: Re: More Inquiries
  • Next by thread: Re: Re: More Inquiries