Re: Re: Integration Program of Burr Distribution

*To*: mathgroup at smc.vnet.net*Subject*: [mg91577] Re: [mg91565] Re: [mg91545] Integration Program of Burr Distribution*From*: Darren Glosemeyer <darreng at wolfram.com>*Date*: Thu, 28 Aug 2008 03:16:24 -0400 (EDT)*References*: <200808260730.DAA23794@smc.vnet.net> <200808271044.GAA21323@smc.vnet.net>

Daniel Lichtblau wrote: > Dewi Anggraini wrote: > >> Dear All >> >> I'd like to calculate the upper specification limit of Burr Distribution from my data below: >> >> 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}; >> >> with the initialisation pdf of Burr as follows: >> >> BurrDistribution[x1_, c_, >> k_] := (c*k)*(x1^(c - 1)/(1 + x1^c)^(k + 1)) >> pdf = BurrDistribution[x1, c, k] >> >> I already got the parameter value of c and k through MLE: >> >> logl = Plus @@ Log[pdf] >> maxlogl = FindMinimum[{-logl, c > 0 && k > 0}, {c, 1}, {k, 2}, >> MaxIterations -> 1000] >> mle = maxlogl[[2]] >> >> However, when I wanna calculate the integration of Burr with the value of c=37.8115, k=0.0135614 and x1 limit = [8, Infinity] to get the probability of data fall outside the upper specification limit given {8} with this program below: >> >> In[2]:= Integrate[pdf[37.81151579009424, 0.01356141249769735, t], >> {t, 8, Infinity}] >> >> Out[2]= Integrate[pdf[37.81151579009424, 0.01356141249769735, t], >> {t, 8, Infinity}] >> >> That gives me the repetation not a real number. >> >> Why is this happen? Do I make mistake somewhere? Because when I did the similar program towards Gamma Distribution, it works well. >> > > > >> I look forward for your reply and assistance. I highly appreciate the suggestions you contribute to my minor thesis. >> >> Thank You. >> Regards, >> Dewi >> > > Among the obvious problems: > > (1) You defined pdf to be a (very large) list, not a function of three > parameters. > > (2) Your integration occurs at In[2]. This means you could not have > (re)defined all the necessary information for pdf. I would surmise you > used Quit[] or otherwise restarted your kernel, and did not redefine things. > > More subtle: > > (3) I doubt your MLE is what you intend. You minimized a sum of logs of > functions of data values (ordinates, or y values). They are in no way > related to the abscissae, or x values. These, I presume, are meant to be > values 1,2,...,69. Also, you will require a sum of discrepancie measures > (absolute values, or maybe squares). So I think your MLE optimization > should involve an expression such as logtot below. > > pdflist = BurrDistribution[Range[n], c, k]; > logtot = Total[Abs[PowerExpand[Log[pdflist]] - Log[x1]]]; > > Or possibly it could just use > tot = Total[Abs[pdflist - x1]]; > > And you could get the optimal values via > > minlogl = > FindMinimum[{logtot, c > 0 && k > 0}, {c, 1}, {k, 2}, > MaxIterations -> 1000]; > mle = minlogl[[2]] > > If you now define > > pdf[t_, c_, k_] := BurrDistribution[t, c, k] > > (note order of arguments) then you can do quadrature as below. Observe > that one can use mle to substitute numeric values for the parameters; > you do not (and should not) cut/paste from the earlier computation. > > NIntegrate[pdf[t, c, k] /. mle, {t, 8, Infinity}] > > > Daniel Lichtblau > Wolfram Research > > Actually, Dewi's mle is fine. To clarify the optimization criterion, the maximum likelihood estimate is the one that maximizes the likelihood, which is the probability of observing the actual data given values of the distributional parameters (c and k in this case). Assuming independent observations, the likelihood is the product of p(xi | c,k) over all i, where p(xi | c,k) is the pdf at xi given c and k. The mle is then obtained by maximizing this product with respect to c and k, or equivalently (since the probabilities are greater than 0) maximizing the sum of the logs of the p(xi| c,k), or minimizing the negative of this sum. Optimizing the sum of logs rather than the product is common and typically puts the objective function on a nicer scale. Darren Glosemeyer Wolfram Research

**References**:**Integration Program of Burr Distribution***From:*"Dewi Anggraini" <dewi_anggraini@student.rmit.edu.au>

**Re: Integration Program of Burr Distribution***From:*Daniel Lichtblau <danl@wolfram.com>