Re: [Rd] qt for df < 1

From: roger koenker <rkoenker_at_uiuc.edu>
Date: Fri 09 Dec 2005 - 17:36:24 GMT

On Dec 9, 2005, at 10:05 AM, Luke Tierney wrote:

> On Thu, 8 Dec 2005, Peter Dalgaard wrote:
>
>> roger koenker <rkoenker@uiuc.edu> writes:
>>
>>> I was experimenting yesterday with a binomial make.link option
>>> for estimating student t binary response models, tentatively
>>> called gossit, and I noticed eventually that the R qt function
>>> doesn't
>>> like df < 1. Vaguely recalling that Splus didn't seem to mind such
>>> weirdness, I checked on our soon to be defunct Splus6.2 and
>>> sure enough, it produced plausible answers instead of R's NA's.
>>> Of course, I have no way of judging the quality of these answers,
>>> but I'm curious about whether someone has already looked into
>>> this can of worms.
>>
>> Well the help page has:
>>
>> For 'qt' only values of at least one are currently supported.
>>
>> and someone must have written that...
>>
>> R does have pt for df < 1, so a temporary fix using uniroot() seems
>> doable.
>>
>>
>
> Something like
>
> qqt<-function(p,df) sign(p-0.5)*sqrt(qf(1-2*pmin(p,1-p),1,df))
>
> seems to do reasonably, at least in terms of consistency with pt, down
> to 0.2 or mayby 0.1 df based on
>
> f<-function(d, n = 101) {
> x<-seq(0, 1, len = n)
> max(abs(pt(qqt(x, d), d) - x))
> }
> plot(function(d) log10(sapply(d, f)), .01,1)
>
> luke

Splus6.2 seems to be using just this approach based on similar testing.

Pure schadenfreude makes it hard to resist mentioning that in Splus6.2 qf(0, df1, df2) gives -Inf rather than 0, which caused some difficulties initially with my
attempt to replicate the comparison.

Roger

>
> --
> Luke Tierney
> Chair, Statistics and Actuarial Science
> Ralph E. Wareham Professor of Mathematical Sciences
> University of Iowa Phone: 319-335-3386
> Department of Statistics and Fax: 319-335-3017
> Actuarial Science
> 241 Schaeffer Hall email: luke@stat.uiowa.edu
> Iowa City, IA 52242 WWW: http://www.stat.uiowa.edu



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat Dec 10 04:41:22 2005

This archive was generated by hypermail 2.1.8 : Fri 09 Dec 2005 - 21:21:08 GMT