[Rd] Problem with fitdistr for gamma in R 2.2.0
Prof Brian Ripley
ripley at stats.ox.ac.uk
Tue Dec 6 12:22:27 CET 2005
The problem lies in dgamma, the function which is stated to be implicated.
See the following NEWS item:
o [dpqr]gamma now returns NaN for an invalid 'shape' parameter
(rather than throw an error), for consistency with other
distribution functions.
That setting lower/upper works is just intelligence in the function to
choose an appropriate method.
On Thu, 17 Nov 2005, P Ehlers wrote:
> I think the problem may lie with fitdistr().
> Specifically, replacing the code in fitdistr.R (VR_7.2-20)
> (line 137 to end) with the code in VR_7.2-8 (line 92 to end)
> seems to handle
>
> fitdistr(otm, "gamma")
>
> just fine. But I haven't done much testing.
>
> Peter Ehlers
>
> Peter Dalgaard wrote:
>> P Ehlers <ehlers at math.ucalgary.ca> writes:
>>
>>
>>> Gregor,
>>>
>>> fitdistr(otm, "gamma", method="L-BFGS-B")
>>>
>>> works for me (on WinXP). Or you could specify "lower = 0".
>>
>>
>> The really odd thing is that it even works with
>>
>>
>>> fitdistr(otm, "gamma",lower=-Inf)
>>
>> shape rate
>> 1.03081094 0.18924370
>> (0.09055117) (0.02117350)
>>
>> or even
>>
>>
>>> fitdistr(otm, "gamma",upper=Inf)
>>
>> shape rate
>> 1.03081094 0.18924370
>> (0.09055117) (0.02117350)
>>
>>
>> Also
>>
>>
>>> fitdistr(otm, "gamma",control=list(parscale=c(.1,.1)))
>>
>> shape rate
>> 1.03079500 0.18923897
>> (0.09055106) (0.02117363)
>>
>> and quite amusingly:
>>
>>
>>
>>> fitdistr(otm, "gamma",method="BFGS",lower=0)
>>
>> shape rate
>> 1.03081096 0.18924371
>> (0.09055118) (0.02117350)
>> Warning message:
>> bounds can only be used with method L-BFGS-B in: optim(x = c(0.059610966029577, 0.0591496321922168, 0.14, 0.18,
>>
>>> fitdistr(otm, "gamma",method="CG",lower=0)
>>
>> shape rate
>> 1.03081096 0.18924371
>> (0.09055118) (0.02117350)
>> Warning message:
>> bounds can only be used with method L-BFGS-B in: optim(x = c(0.059610966029577, 0.0591496321922168, 0.14, 0.18,
>>
>> whereas the same calls without the dysfunctional lower= gives the
>> warning about `shape` needing to be positive.
>>
>> This probably all indicates that something inside optim() is broken.
>>
>>
>>
>>
>>> I no longer have 2.1.0 running, so I don't know why this
>>> wasn't needed in 2.1.0.
>>>
>>> "R version 2.2.0, 2005-10-24"
>>> MASS version: 7.2-20
>>>
>>> -peter
>>>
>>> Gorjanc Gregor wrote:
>>>
>>>> Dear R developers,
>>>>
>>>> I have encountered strange behaviour of fitdistr for gamma in recent R
>>>> build i.e. 2.2.0. I have attached the code for data at the end of this mail
>>>> so you can reproduce the problem. In short, I am able to run fitdistr under
>>>> 2.1.0 without problems, while I get the following error under 2.2.0
>>>> (Version 2.2.0 Patched (2005-11-15 r36348))
>>>>
>>>>
>>>>
>>>>> fitdistr(otm, "gamma")
>>>>
>>>> Error in densfun(x, parm[1], parm[2], ...) :
>>>> 'shape' must be strictly positive
>>>>
>>>> The results on 2.1.1 (Version 2.1.1 (2005-06-20)) are
>>>>
>>>>
>>>>
>>>>> fitdistr(otm, "gamma")
>>>>
>>>> shape rate
>>>> 1.030667 0.189177
>>>> (0.090537) (0.021166)
>>>>
>>>> Platform: Windows XP
>>>>
>>>> Thank you in advance for your effort on this remarkable tool!
>>>>
>>>> Here is the data for above problem/results:
>>>>
>>>> "otm" <-
>>>> c(0.059610966029577, 0.0591496321922168, 0.14, 0.18, 0.24, 0.25,
>>>> 0.270071982912719, 0.270758049933706, 0.269911804412492, 0.280138451903593,
>>>> 0.279787947586738, 0.279429937571753, 0.3, 0.320746235495899,
>>>> 0.319553311037365, 0.51, 0.54, 0.56, 0.6, 0.609812622915953,
>>>> 0.609198293855879, 0.64, 0.69, 0.74, 0.76, 0.770972826186568,
>>>> 0.769288654833566, 0.78, 0.789181584270671, 0.78991363293305,
>>>> 0.8, 0.89, 0.900691718998831, 0.8991656800583, 0.92, 0.93, 0.94,
>>>> 1.01, 1.02, 1.13, 1.18, 1.26, 1.29, 1.33, 1.42, 1.43, 1.47, 1.47940529614314,
>>>> 1.47920716832764, 1.6, 1.61, 1.63, 1.68938231960637, 1.6894849291523,
>>>> 1.82, 1.88088044053270, 1.8792804789003, 1.89, 1.92, 2, 2.04,
>>>> 2.07, 2.12, 2.17, 2.18, 2.22, 2.23, 2.27, 2.28, 2.3, 2.32092240267433,
>>>> 2.31912300181622, 2.38, 2.39, 2.43, 2.46, 2.51, 2.52, 2.55, 2.56,
>>>> 2.61, 2.66091404781397, 2.6595832825806, 2.67, 2.7, 2.77, 2.8,
>>>> 2.81, 2.86, 2.87, 2.93, 3.01, 3.05, 3.14, 3.15, 3.17, 3.18, 3.24,
>>>> 3.26, 3.33, 3.44, 3.45, 3.52, 3.55, 3.63, 3.73, 3.9, 4, 4.01,
>>>> 4.04, 4.13, 4.15934497380769, 4.16094719917513, 4.3, 4.33, 4.34,
>>>> 4.66, 4.76, 4.82, 4.83, 4.89, 4.92, 5.06, 5.14, 5.16, 5.26, 5.31,
>>>> 5.36, 5.48, 5.66, 5.79, 5.8, 5.85, 5.87, 5.92952534468565, 5.92962284128508,
>>>> 6.04, 6.11, 6.13, 6.16, 6.19, 6.42, 6.66, 6.69, 7.11, 7.16, 7.29,
>>>> 7.3, 7.31, 7.33, 7.72, 7.82, 7.87, 7.91, 8.01, 8.17, 8.45, 8.49,
>>>> 8.73, 8.86, 8.95, 9, 9.05, 9.13, 9.22, 9.52, 9.82, 9.88, 9.91,
>>>> 9.99, 10.03, 10.4, 10.59, 10.83, 11.06, 11.64, 11.85, 12.02,
>>>> 12.4, 12.64, 12.96, 13.44, 14.06, 14.07, 14.37, 15.4, 15.6, 15.92,
>>>> 16.23, 16.6, 16.97, 17.06, 17.8, 18.69, 18.73, 19.2, 19.51, 19.54,
>>>> 20.57, 21.05, 22.23, 27.02)
>>>>
>>>> Lep pozdrav / With regards,
>>>> Gregor Gorjanc
>>>>
>>>> ----------------------------------------------------------------------
>>>> University of Ljubljana
>>>> Biotechnical Faculty URI: http://www.bfro.uni-lj.si/MR/ggorjan
>>>> Zootechnical Department mail: gregor.gorjanc <at> bfro.uni-lj.si
>>>> Groblje 3 tel: +386 (0)1 72 17 861
>>>> SI-1230 Domzale fax: +386 (0)1 72 17 888
>>>> Slovenia, Europe
>>>> ----------------------------------------------------------------------
>>>> "One must learn by doing the thing; for though you think you know it,
>>>> you have no certainty until you try." Sophocles ~ 450 B.C.
>>>>
>>>> ______________________________________________
>>>> R-devel at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>
>>> Peter Ehlers
>>> University of Calgary
>>>
>>> ______________________________________________
>>> R-devel at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>
>>
>>
>
> --
> Peter Ehlers
> Department of Mathematics and Statistics
> University of Calgary, 2500 University Dr. NW
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-devel
mailing list