Simple addition to random module - Student's t

Thomas Philips tkpmep at gmail.com
Wed Sep 2 13:15:57 EDT 2009


On Sep 2, 1:03 pm, Thomas Philips <tkp... at gmail.com> wrote:
> On Sep 2, 12:28 pm, Mark Dickinson <dicki... at gmail.com> wrote:
>
> > On Sep 2, 2:51 pm, Thomas Philips <tkp... at gmail.com> wrote:
>
> > > def student_t(df):         # df is the number of degrees of freedom
> > >     if df < 2  or int(df) != df:
> > >        raise ValueError, 'student_tvariate: df must be a integer > 1'
>
> > By the way, why do you exclude the possibility df=1 here?
>
> > --
> > Mark
>
> I exclude df=1 hereBecause the variance is then infinite (in fact, the
> distribution is then Cauchy). That said, your point is well taken;
> allowing df=1 makes the Cauchy distribution available to users of
> random, in much the same way as the Gamma makes the Chi-squared
> available. Here's the revised code:
>
> def student_tvariate(df):         # df is the number of degrees of
> freedom
>     if df < 1  or int(df) != df:
>         raise ValueError, 'student_tvariate: df must be a positive
> integer'
>
>     x = random.gauss(0, 1)
>     y = random.gammavariate(df/2.0, 2)
>
>     return x / (math.sqrt(y/df))
>
> I'll follow your suggestion, add in documentation and submit it to
> bugs.python.com. Thanks for your guidance.
>
> Thomas Philips

Mark,

I mis-spoke - the variance is infinite when df=2 (the variance is df/
(df-2), and you get the Cauchy when df=2. So I'm going to eat humble
pie and go back to

def student_tvariate(df):         # df is the number of degrees of
freedom
    if df < 2  or int(df) != df:
       raise ValueError, 'student_tvariate: df must be a integer > 1'

    x = random.gauss(0, 1)
    y = random.gammavariate(df/2.0, 2)

    return x / (math.sqrt(y/df))



I made the mistake because the denominator is  equivalent to the
square root of the sample variance of df normal observations, which in
turn has df-1 degrees of freedom...............Oh well,  I think it's
much easier to apologize than to rationalize my error!!!!


Sincerely

Tom



More information about the Python-list mailing list