Re: [ferret_users] Statistics / student-t distribution

[Andrew Wittenberg <andrew.wittenberg@xxxxxxxx>]

Ansley et al.,

The external function STUDENT_T_CUTOFF could do with more
documentation.  To clarify, it doesn't return a quantile of the
t-distribution with nf degrees of freedom.  Instead,
STUDENT_T_CUTOFF(P,nf) returns the P/100 critical value (upper limit)
for the two-sided, one-sample t-test that the mean of a sample of nf
independent normally-distributed observations differs from zero.  Note
that STUDENT_T_CUTOFF is rather limited: it only accepts P-values of
90%, 95%, or 99%.

Given a quantile q and degrees of freedom df, the mapping onto
STUDENT_T_CUTOFF is:

   P = 100*(1-2*(1-q))
   nf = df+1

and the reverse mapping is

   q = 1-(1-P/100)/2
   df = nf-1

See e.g. "One-sample t-test" under

   http://en.wikipedia.org/wiki/Student%27s_t-test

To obtain an arbitrary t-quantile, I'd recommend calling R or Python
(scipy.stats) for now, until PyFerret arrives with its powerful
statistics functions.  With a slight modification to Paulo's
suggestion, we can remove the dependence on "littler" ("little R") and
just call R directly as a one-liner.  E.g. to compute the 95th
percentile of the t-distribution with 10 degrees of freedom:

   let qt = STRFLOAT(SPAWN("echo 'cat(qt(0.95, df=10))' | R --slave --vanilla"))

Andrew

PS. For those interested, here's a comparison of STUDENT_T_CUTOFF with
the t-quantile (qt) function from R:

        NOAA/PMEL TMAP
        FERRET v6.72
        Linux 2.6.18-274.el5PAE 32-bit - 09/13/11
        21-Mar-12 14:02

yes? let q = .95; let df = 10
yes? list student_t_cutoff(100*(1-2*(1-q)),df+1), strfloat(spawn("echo
'cat(qt(`q`, df=`df`))' | R --slave --vanilla"))
             X: 1
 Column  1: STUDENT_T_CUTOFF(100*(1-2*(1-Q)),DF+1)
 Column  2: STRFLOAT(SPAWN("echo 'cat(qt(0.95, df=10))' | R --slave --vanilla"))
       (C001,V00  (C001,V008)
I / *:      1.812   1.812

On Wed, Mar 21, 2012 at 12:21 PM, Ansley Manke <ansley.b.manke@xxxxxxxx> wrote:
> Hi all,
> Ferret does include an external function, student_t_cutoff, which may
> provide some of what you're looking for:
>
> yes?  sh func student*
> STUDENT_T_CUTOFF(P,nf)
>    Return student-t cutoff
>    P: Confidence Limit
>    nf: nf (Number in sample)
>
> Other statistics capabilities would be straitforward to write as external
> functions.  Our current development, PyFerret, builds Ferret as a Python
> module. This will in effectHi Ansley,

The external function STUDENT_T_CUTOFF could use more documentation.
To clarify, it doesn't return a quantile of the t-distribution with nf
degrees of freedom.  Instead, STUDENT_T_CUTOFF(P,nf) returns the P/100
critical value (upper limit) for the two-sided, one-sample t-test that
the mean of a sample of nf independent normally-distributed
observations differs from zero.  Note that STUDENT_T_CUTOFF is rather
limited, as it only accepts P-values of 90, 95, or 99.

Given a quantile q and degrees of freedom df, the mapping onto
STUDENT_T_CUTOFF is:

   P = 100*(1-2*(1-q))
   nf = df+1

and the reverse mapping is

   q = 1-(1-P/100)/2
   df = nf-1

See e.g. "One-sample t-test" under

   http://en.wikipedia.org/wiki/Student%27s_t-test

To obtain a general t-quantile, I'd recommend using Python
(scipy.stats) or R for now, until PyFerret with its powerful
statistics functions) becomes available.  With a slight modification
to Paulo's suggestion, we can elminate the dependence on "littler"
("little R") and just call R directly.  E.g. to compute the 95th
percentile of the t-distribution with 10 degrees of freedom:

   let qt = STRFLOAT(SPAWN("echo 'cat(qt(0.95, df=10))' | R --slave --vanilla"))

Andrew

PS. For those who are interested, here's a comparison of
STUDENT_T_CUTOFF with the t-quantiles:

        NOAA/PMEL TMAP
        FERRET v6.72
        Linux 2.6.18-274.el5PAE 32-bit - 09/13/11
        21-Mar-12 14:02

yes? let q = .95; let df = 10
yes? list student_t_cutoff(100*(1-2*(1-q)),df+1), strfloat(spawn("echo
'cat(qt(`q`, df=`df`))' | R --slave --vanilla"))
             X: 1
 Column  1: STUDENT_T_CUTOFF(100*(1-2*(1-Q)),DF+1)
 Column  2: STRFLOAT(SPAWN("echo 'cat(qt(0.95, df=10))' | R --slave --vanilla"))
       (C001,V00  (C001,V008)
I / *:      1.812   1.812

On Wed, Mar 21, 2012 at 12:21 PM, Ansley Manke <ansley.b.manke@xxxxxxxx> wrote:
> Hi all,
> Ferret does include an external function, student_t_cutoff, which may
> provide some of what you're looking for:
>
> yes?  sh func student*
> STUDENT_T_CUTOFF(P,nf)
>    Return student-t cutoff
>    P: Confidence Limit
>    nf: nf (Number in sample)
>
> Other statistics capabilities would be straitforward to write as external
> functions.  Our current development, PyFerret, builds Ferret as a Python
> module. This will in effect make mathematical and statistical functions from
> Python libraries available as Ferret functions.
>
> Ansley
>
> On 3/21/2012 8:31 AM, Paulo B. Oliveira wrote:
>>
>> One way is to use a statistical package to do the job. I ended up
>> learning some R (http://www.r-project.org/) basics (the choice was
>> driven by its popularity among my colleagues).
>> When combined with 'littler' ( http://code.google.com/p/littler/ ) you
>> can use spawn in ferret to get the simple statistics.
>>
>> For example, to get student-t 0.95 quantile with df-degrees of freedmon
>> from ferret:
>> yes? def sym df 100
>> yes? let qt95 = `spawn("r -e 'cat(qt((0.95),df=($df)))'")`
>>  !->  DEFINE VARIABLE qt95 = 1.66023
>>
>> If a long data series has to be passed on to R, it will be necessary to
>> write an R script, eventually using the RNetCDF
>> ( http://cran.r-project.org/web/packages/RNetCDF/index.html )
>>  ncdf ( http://cran.r-project.org/web/packages/ncdf/index.html ) or
>> ncdf4 (http://cran.r-project.org/web/packages/ncdf4/index.html )
>> packages.
>>
>>
>> Regards,
>>
>> Paulo
>> ----------------------------------------------------
>> On Mon, 2012-03-19 at 09:02 +0100, Hein Zelle wrote:
>>>
>>> Neil Swart wrote:
>>>
>>>> I wondering how other Ferreters compute their statistics, and if
>>>> anyone has a TINV script or similar already written?
>>>
>>> I have not done this yet, but for an upcoming job I'll need to
>>> determine similar bounds (quantiles) from an ensemble distribution.
>>> I can do this based on the std deviation, assuming a gaussian
>>> distribution, but that doesn't always hold.
>>>
>>> I'll be very interested to hear if it's possible to do this using
>>> external functions or with different tricks.
>>>
>>> Kind regards,
>>>      Hein Zelle
>>>
>>>
>>
>

 make mathematical and statistical functions from
> Python libraries available as Ferret functions.
>
> Ansley
>
> On 3/21/2012 8:31 AM, Paulo B. Oliveira wrote:
>>
>> One way is to use a statistical package to do the job. I ended up
>> learning some R (http://www.r-project.org/) basics (the choice was
>> driven by its popularity among my colleagues).
>> When combined with 'littler' ( http://code.google.com/p/littler/ ) you
>> can use spawn in ferret to get the simple statistics.
>>
>> For example, to get student-t 0.95 quantile with df-degrees of freedmon
>> from ferret:
>> yes? def sym df 100
>> yes? let qt95 = `spawn("r -e 'cat(qt((0.95),df=($df)))'")`
>>  !->  DEFINE VARIABLE qt95 = 1.66023
>>
>> If a long data series has to be passed on to R, it will be necessary to
>> write an R script, eventually using the RNetCDF
>> ( http://cran.r-project.org/web/packages/RNetCDF/index.html )
>>  ncdf ( http://cran.r-project.org/web/packages/ncdf/index.html ) or
>> ncdf4 (http://cran.r-project.org/web/packages/ncdf4/index.html )
>> packages.
>>
>>
>> Regards,
>>
>> Paulo
>> ----------------------------------------------------
>> On Mon, 2012-03-19 at 09:02 +0100, Hein Zelle wrote:
>>>
>>> Neil Swart wrote:
>>>
>>>> I wondering how other Ferreters compute their statistics, and if
>>>> anyone has a TINV script or similar already written?
>>>
>>> I have not done this yet, but for an upcoming job I'll need to
>>> determine similar bounds (quantiles) from an ensemble distribution.
>>> I can do this based on the std deviation, assuming a gaussian
>>> distribution, but that doesn't always hold.
>>>
>>> I'll be very interested to hear if it's possible to do this using
>>> external functions or with different tricks.
>>>
>>> Kind regards,
>>>      Hein Zelle
>>>
>>>
>>
>