WOLFRAM

Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2008

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

Search the Archive

Re: Re: Integration Program of Burr Distribution

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
  • Prev by Date: Re: problem with using external functions
  • Next by Date: Re: efficiently adding many 2D Gaussians
  • Previous by thread: Re: Integration Program of Burr Distribution
  • Next by thread: Catching messages