Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

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

Search the Archive

Re: Please help me.....

  • To: mathgroup at smc.vnet.net
  • Subject: [mg91101] Re: [mg91072] Please help me.....
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Wed, 6 Aug 2008 05:05:22 -0400 (EDT)
  • References: <200808050803.EAA09872@smc.vnet.net>

Dewi Anggraini wrote:
> Hi all
> 
> I'm an international studentat RMIT University Melbourne. Currently, I'm 
> doing my master degree in Statistics and Operations Research. This
> following 6 months will be my last semester to finish my master degree.
> Therefore, I have been doing minor thesis about "estimating unknown
> parameter of Burr distribution" by using Mathematica.
> 
> I was wondering if this forum could assist me to detect where I got wrong
> with my program since it works for Burr distribution in some data but not
> with my data. Additionally, this program (I attached my program along with
> my data) also works for non-normal distribution, such as Gamma and Weibull
> since my data are non-normally distributed and they are closely to gamma
> and weibull distribution.
> 
> The data are treatment time of cervical cancer patients in a hospital.
> Frequently, they come to the hospital in the late satge of cancer, thus
> the minimum time for the treatment is 2 days and the maximum treatment on
> the data is 21 days.
> 
> The following is my program to run MLE of Burr distribution. However,
> it comes up with "comment" results in finding the coefficient of unknown
> parameters when I run it.
> 
> 
> n = 69;
> 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};
> BurrDistribution[x2_, c_,
>   k_] := (c*k)*(x2^(c - 1)/(1 + x2^c)^(k + 1))
> pdf = BurrDistribution[x2, c, k]
> logl = Plus @@ Log[pdf]
> maxlogl = FindMinimum[-logl, {c, 1}, {k, 2}]
> mle = maxlogl[[2]]
> 
> 
> 
> Please assist me in finding the problem I face now. This is very important
> for my thesis.
> 
> I look forward for any comments and advices, thank you.
> 
> Regards,
> Dewi

I am not sure what is the issue but I suspect you refer to the warning 
messages. They indicate that some log arguments probably became 
negative, hence caused complex values to appear. There are various ways 
to address that and make the job easier.

One is to expand the logs, under the assumption that the parameters c 
and k are nonnegative. This can be done as below.

logl = PowerExpand[Total[Log[pdf]], Assumptions->{c>=0,k>=0}];

Also you can restrict FindMinimum by giving constraints, or instead use 
NMinimize and likewise give constraints. Either of the calls below will 
do this.

In[23]:= maxlogl = FindMinimum[{-logl,c>=0,k>=0}, {c,1}, {k,2}]

FindMinimum::eit:
                                                                  -6
    The algorithm does not converge to the tolerance of 4.80622 10   in 500
      iterations. The best estimated solution, with feasibility 
residual, KKT
                                                       -12
      residual or complementary residual of {7.10845 10   , 0.0000304113,
                -12
      5.93765 10   }, is returned.

Out[23]= {226.436, {c -> 20.2248, k -> 0.0285432}}

In[24]:= maxloglb = NMinimize[{-logl,c>=0,k>=0}, {{c,1,2}, {k,2,4}}]

Out[24]= {226.436, {c -> 15.1598, k -> 0.0380797}}

You might observe that they give the same minimum, to six digits, but 
give quite different values for the parameters {c,k}. Whether this is 
reasonable is something I cannot address.

Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: Derivative of Dot[]
  • Next by Date: RE: Re: Workaround for an unexpected behavior of Sum
  • Previous by thread: Re: Please help me.....
  • Next by thread: Re: Please help me.....