slow running betainc function

From: Will Hobbs <whobbs_at_nyahnyahspammersnyahnyah>
Date: Wed, 01 Aug 2007 15:09:55 PDT

('binary' encoding is not supported, stored as-is) Hi

I want to do a paired-difference test between two model simulations, and
have a function (below) that will perform and one-sample t-test (pasted)
below. With the 'tval' flag set to False so that the function returns the
2-tailed probability, the probility calculation using the incomplete beta
function is painfully slow. The same is true if I use the new 'student_t'
function as well. Am I doing something wrong? Is there anything I can do
to make this more efficient?

For the record, the input array sizes are (diff= 3, lat= 25, lon= 128)

Many thanks

Will

;******************************************************************************
;ttest_1samp
;***************************************************************

;Performs one-sample Student's t-test, typically for paired-difference test

;prototype
  function ttest_1samp(x:numeric, var:numeric, N:integer, c[1]:numeric,
tval[1]:logical)

;description

;x : array of any dimension of means
;var : array of variances, must be same dimensions as x
;N : scalar or array of sam size as x of sample size
;c : scalar value against which x is compared. Typically 0.
;tval : if True, t values are returned, if false probabilities are returned.

;Null hypothesis is that x - c = 0

  local xdim, vdim, ndim, t, df, prob

  begin

        ;test that x, var are same dimensions

        xdim = dimsizes(x)
        vdim = dimsizes(var)

        if (any(xdim.ne.vdim)) then
           print("fatal: ttest_1samp. Mean and variance arrays must have same
dimensions")
           exit
        end if

        ;test that N is scalar or correctly sized array

        ndim = dimsizes(N)
        

        if (sum(ndim).gt.1.and.any(ndim.ne.vdim)) then
                print("fatal: ttest_1samp. Sample size must be a scalar, or have same
dimensions as x")
                exit
        end if
                 

;calculate t-value

        t = (x - c) / (( var/N)^0.5 )

        if (tval.eq.True) then

                return(t)

        else
                df = N-2

;this is the (really) slow bit!
                prob = betainc( df/(df+t^2), df/2, conform(t, 0.5, -1))
                return(prob)
                
        end if
        
end

============================
Will Hobbs
PhD Candidate
Room A113
UCLA Department of Geography
1255 Bunche Hall
Los Angeles CA 90095-1524
Tel: 310 663 2631

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Aug 01 2007 - 16:09:55 MDT