[R] Strange behavior using .Fortran

[Prof Brian Ripley]

You need to read the assumptions:

      SUBROUTINE EXCHNG (X, M, Y, N, SX, SY)
C
C        ALGORITHM AS 304.2 APPL.STATIST. (1996), VOL.45, NO.3
C
C        Exchanges the sample data.  Assumes both X and Y have been
C        previously dimensioned to at least max(M, N) elements

and you have not done this!  Try

fisher.ts.test <- function(x,y)
{
  nx <- length(x)
  ny <- length(y)
  xx <- yy <- double(max(nx, ny))
  xx[1:nx] <- x
  yy[1:ny] <- y
  .Fortran("fisher", xx, as.integer(nx),
           yy, as.integer(ny),
           total=integer(1),
           possib=integer(1),
           P=double(1),
           ierr=integer(1)
           )
}

I also suggest you compile by 

R SHLIB fisher.f

as you will then get, correctly, a shared library with an .so extension.


> From: "Uli Flenker; Raum 704" <uli at biochem.dshs-koeln.de>
> To: R-Mailing-List <r-help at stat.math.ethz.ch>
> Subject: [R] Strange behavior using .Fortran
> 
> Hello everybody,
> 
> I observed some strange behavior of R (0.64.1/Linux 2.0.25) when trying
> to use FORTRAN-libraries.
> 
> I downloaded "ALGORITHM AS 304" from StatLib. All I changed  at
> the code, was to substitute "REAL"-routines by "DOUBLE PRECISION".
> 
> Compiling went o.k. ...
> 
> > g77 -fpic -O2 -shared -o rtest.o fisher.f
> 
>  ("g77 --version" gives 2.7.2.1)
> 
> 
> Here are the first lines of the code from AS 304, "FISHER" being the
> routine to be used from within R:
> 
>       SUBROUTINE FISHER (X, M, Y, N, TOTAL, POSSIB, P, IFAULT)
> C
> C        ALGORITHM AS 304.1 APPL.STATIST. (1996), VOL.45, NO.3
> C
> C        Fisher's non-parametric randomization test for two small
> C        independent random samples
> C
>       INTEGER M, N, TOTAL, POSSIB, IFAULT
>       DOUBLE PRECISION X(*), Y(*), P
> 
> C
> 	... etc., lots of code ...
> 
> 
> 
> Here is the prototype of R-code that I used ...
> 
> > dyn.load("/path-to-lib/rtest.o")
> 
> > fisher.ts.test<-function(x,y){
> 
>    nx<-length(x)
>    ny<-length(y)
> 
>          
>    ftr.out<-..Fortran("fisher",
>                      as.real(x),
>                      as.integer(nx),
>                      as.real(y),
>                      as.integer(ny),
>                      as.integer(0),
>                      as.integer(0),
>                      as.real(0.0),
>                      as.integer(0))
> 
> return(ftr.out)
> }
> 
> 
> Loading the code was no problem. Now, here is what I obtained when using
> it:
> 
> > a<-c(18,19,25)
> > b<-c(20,29,29,30,31,31,31)
> 
> > fisher.ts.test(a,b)
> [[1]]
> [1] 18 19 25
> 
> [[2]]
> [1] 3
> 
> [[3]]
> [1] 20 29 29 30 31 31 31
> 
> [[4]]
> [1] 7
> 
> [[5]]
> [1] 1
> 
> [[6]]
> [1] 120
> 
> [[7]]
> [1] 0.01666667
> 
> [[8]]
> [1] 0
> 
> ... which are the correct results.
> 
> 
> But invoking the thing exchanging "a" and "b" gives the following strange
> things ...
> 
> > fisher.ts.test(b,a)
> [[1]]
> [1] 31 31 31 30 31 31 31
> 
> [[2]]
> [1] 3
> 
> [[3]]
> [1] 20 29 29
> 
> [[4]]
> [1] 0
> 
> [[5]]
> [1] 0
> 
> [[6]]
> [1] 1
> 
> [[7]]
> [1] 1
> 
> [[8]]
> [1] 0
> 
> 
> Furthermore ...
> 
> > x<-round(rnorm(3,10),2)
> > y<-round(rnorm(7,5),2)
> > fisher.ts.test(x,y)
> Process R segmentation fault at Thu Jul 15 19:08:01 1999
> 
> 
> This could not be observed when using a similar routine written in
> Xlisp-Stat:
> 
> (dyn-load "/path-to-lib/rtest.o")
> 
> (defun tse-test (x y)
>   (let* (
> 	 (out-list (call-fsub "fisher" (float x) (length x)
> 			               (float y) (length y)
> 				        0 0 (float 0) 0))
> 	 
> 	 (pval   (elt out-list 6))
> 	 (poss   (elt out-list 5))
> 	 (par    (elt out-list 4))
> 	 )
> 
> 
> (format T "\nFISHER's non-parametric Randomization Test
> for two small independent Samples
> 
> Z             :~8,d          
> Possibilities :~8,d
> p (one-sided) :~8,5,f\n" (elt par  0)(elt poss 0)(elt pval 0)) 
> )
> )
> 
> This function worked fine under any circumstances until now.
> (Xlisp-Stat 3.52.9)
> 
> A small FORTRAN-Program that allowed to type in the data
> interactively also didn't care about the order of the size
> or sum of the samples.
> 
> All this makes me conclude the reason is somewhere in R.

No, the problem is in your reading of the conditions of the Fortran code.
You were just lucky that this works in Xlisp-Stat, as you were
writing off the end of a Fortran array.

-- 
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 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._