help links subscribe changes refresh edit

	SandBox
SandBoxFileUploadTest ←	↑	→ SandBoxFormalProduct

SandBoxFisher

Edit detail for SandBoxFisher revision 1 of 11

	1 2 3 4 5 6 7 8 9 10 11
Editor: Mark Clements Time: 2009/03/16 01:32:09 GMT-7
Note:

changed:
-
How useful are the different CAS languages for implementing numerical routines? Prompted by a comparison of R and C for implementing Fisher's exact test for 2x2 tables (http://fluff.info/blog/arch/00000172.htm), I thought that it would be interesting to implement this particular test in Spad, Boot and Common Lisp (see below). Each set of code was required to implement a univariate root finder and the hypergeometric distribution to calculate the p-value under different alternatives, together with the 95% confidence interval and the maximum likelihood estimator for the odds ratio. The reference implementation is R, where the code and output would be:

\begin{verbatim}
> fisher.test(matrix(c(10,10,10,20),nrow=2))

        Fisher's Exact Test for Count Data

data:  matrix(c(10, 10, 10, 20), nrow = 2) 
p-value = 0.2575
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
 0.5383996 7.4363242 
sample estimates:
odds ratio 
  1.971640 
\end{verbatim}

To summarise, all three languages (Spad, Boot and Common Lisp) provide arbitrary length integers and fractions, ensuring that the hypergeometric distribution was straightforward to implement. The Lisp and R implementations were very similar, which is not surprising, given that the two languages are closely related. In contrast, the nested functions in Spad seemed clunky, with the requirement to use the #1 and #2 argument references, although I did appreciate Spad's lexical scoping and facility to fall back to an algebraic analysis. An Aldor implementation would be cleaner than the Spad implementation. Boot's implementation was initially difficult, as I was unclear how to pass values to the nested functions with using function argument (which was not possible for the univariate root finder :-(). The use of the "$" prefix on derived variables seemed clunky - and I was surprised to find that function arguments were lexically scoped. Finally, my version of Fricas:Boot defaulted to single precision floats, which caused problems with precision. Type specification in OpenAxiom:Boot and more recent versions of Fricas would negate the need for explicit coercion to double-floats.

This begs the question: when would one use any of these languages for mixed numerical/algebraic analysis? First, Boot is the least likely to be used, although it does play closely and well with Common Lisp (a la Maxima). One could code for numerical analysis in Boot and Common Lisp - however Boot's lack of lexical scoping may be a detraction. Moreover, by my understanding, Lisp or Boot are unable to evaluate Spad or Axiom functions. Second, for an R user, Spad's type system seems fussy and debugging can be slow. The lack of ability to call Spad functions from Lisp or Boot is a major restriction, requiring that a large body of code in Lisp (or Fortran via f2cl), such as for optimisation, would need to be hand translated to Spad (see [SandBoxMLE] for an example). As a final note, Maxima would seem to address many of these concerns - although it lacks an elegant type system and I haven't learnt to appreciate its scoping rules :-).


\begin{spad}
)abbrev package TESTP TestPackage
R ==> Float
I ==> Integer
fisherRec ==> Record(PValue:R, CI:List R, Estimate:R)
TestPackage: with
   ridder: (R->R,R,R) -> R
   msign: (R,R) -> R
   choose:(I, I)->Fraction I 
   --chooseNew:(Integer, Integer)->Fraction Integer 
   dhyper:(I, I, I, I)->Fraction Integer
   phyper:(I, I, I, I, Boolean)->Fraction Integer
   fisherTest:(I,I,I,I, String, R, Boolean, R)-> fisherRec
   testTolerance:(R, R, R)->Boolean   
   test1: () -> Boolean
   test2: () -> Boolean
   test3: () -> Boolean
   test4: () -> Boolean
   test5: () -> Boolean
   test6: () -> Boolean
   test7: () -> Boolean
   test8: () -> Boolean
   test9: () -> Boolean
   test10: () -> Boolean
   alltests: () -> List Boolean
  == add
   import TrigonometricFunctionCategory -- for test1()
   --import OrderedCompletion(Float) for plusInfinity()$OrderedCompletion(Float)
   ridder(func, x1, x2) ==
    eps:= 1.0e-16::R
    maxit:= 30::Integer
    --verbose:= false
    fl:R := func x1
    fh:R := func x2
    xl:R := x1
    xh:R := x2
    ans:R := -1.11e30::R
    xnew:R := 0.0e0::R
    iterNum:= 0::Integer
    if fl=0.0::R then return x1
    else if fh=0.0::R then return x2
    else if (fl*fh) > 0.0::R then error "Initial points are not either side of zero."
    --if (fl*fh) < 0.0 then
    else repeat
		xm:= 0.5::R *(xl+xh)
		fm:= func xm
		ss:= sqrt((fm*fm) - (fl*fh))
		if ss =0.0::R then return ans
		xnew:= xm + (((xm - xl) * (if (fl>fh) then 1.0::R else -1.0::R) * fm) / ss)
		if abs(xnew-ans) <= eps then return ans
		ans:= xnew 
		fnew:= func ans
		if fnew=0.0::R then return ans
		if msign(fm,fnew) ~= fm then
		    xl:= xm 
		    fl:= fm 
		    xh:= ans 
		    fh:= fnew
		else if msign(fl, fnew) ~= fl then
		    xh:= ans 
		    fh:= fnew
		else if msign(fh, fnew) ~= fh then
		    xl:= ans 
		    fl:= fnew
		iterNum:=iterNum+1::Integer
		if iterNum >=maxit then 
			error "Maximum iterations exceeded"
		--if verbose then FORMAT(true,"~,8f ~,8f ~,8f ~,8f~%", xl, xh, fl, fh)$Lisp
		if abs(xh-xl) <= eps then return ans
   msign(x, y) ==
    (abs x) * (if y>0.0::R then 1.0::R else if y<0.0::R then -1.0::R else 0.0::R)
   choose(n, x) ==
    total:Fraction Integer := 1/1
    for denom in 1..x repeat
	total:=total*((n-denom+1)/(denom))::Fraction Integer
    return total
   --chooseNew(n, x) == product((n-i+1)::Fraction Integer/i::Fraction Integer,i=1..x)
   dhyper(x, m, n, k) ==
    choose(m, x) * choose(n, k - x) / choose(m + n, k)
   phyper(x, m, n, k, lowerTail) ==
    i:PositiveInteger
    --total:Fraction Integer:=0/1
    if lowerTail then 
	reduce("+",[dhyper(i, m, n, k) for i in 1..x])
    else 
	reduce("+",[dhyper(i, m, n, k) for i in (x+1)..k])
   fisherTest(a,b,c,d, alternative, OR, confInt, confLevel) ==
	m:I := a+c -- first column
	n:I := b+d -- second column
	k:I := a+b -- first row
	x00:I := a
	lo:I := max(0, k-n)
	hi:I := min(k, m)
	support:List I := [i for i in lo..hi]
	logdc:List R:= [log(dhyper(i, m, n, k)::R) for i in support]
	doubleEps:R := 1.0e-50::R 
	plusInfinity:R := 1.0e6::R -- arbitrary
        dnhyper:(R->List R) :=
	     ncp:R := #1
	     d:List R := [logdc(i)+log(ncp)*support(i)::R for i in 1..#logdc]
	     maxd:R := reduce(max,d)
	     d2:List R :=[exp(di-maxd) for di in d]
	     sumd2:R := reduce("+",d2)
	     [d2i/sumd2 for d2i in d2]
        mnhyper:(R->R) :=
             ncp:R := #1
	     if ncp=0.0::R then lo::R
 	     --else if ncp=%plusInfinity then hi::R
	     else 
		d:List R := dnhyper(ncp)
 		reduce("+",[support(i)::R*d(i) for i in 1..#d])
	pnhyper:((Integer,R,Boolean)->R) :=
	     q:I := #1
	     ncp:R := #2
	     upperTail:Boolean := #3
 	     if ncp=1.0 then 
 		if upperTail then phyper(q-1, m, n, k, false)::R
		else phyper(q, m, n, k, true)::R
 	     else if ncp=0.0 then
 		if upperTail then 
 		    if q<=lo then 1.0::R else 0.0::R
 		else if q>=lo then 1.0::R else 0.0::R
-- 	     else if ncp=%plusInfinity then
-- 		if upperTail then 
-- 		    if q<=hi then 1.0::R else 0.0::R
-- 		else if q>= hi then 1.0::R else 0.0::R
 	     else 
		d:List R := dnhyper(ncp)
 		if upperTail then
 		      reduce("+",[d(i) for i in 1..#d | support(i)>=q])
 		else reduce("+",[d(i) for i in 1..#d | support(i)<=q])
	mle:(I->R) :=
	     x:I := #1
	     if x=lo then 0.0::R
	     else if x=hi then plusInfinity
	     else
		mu:R := mnhyper(1.0::R)
		if mu>x::R then 
			f:(R->R) := mnhyper(#1) - x::R
			ridder(f,0,1)
		else if mu<x::R then
			f:(R->R) := mnhyper(1/#1) - x::R
			1/ridder(f,doubleEps,1.0::R)
		else 1.0::R
	ncpU:(I,R)->R :=
	     x:I := #1
	     alpha:R := #2
	     if x=hi then plusInfinity
	     else
		     p:R := pnhyper(x, 1.0::R, false)
		     if p<alpha then 
			f:(R->R) := pnhyper(x,#1,false) - alpha
			ridder(f, 0.0::R, 1.0::R)
		     else if p>alpha then 
				f:(R->R) := pnhyper(x,1/#1,false) - alpha
				1/ridder(f, doubleEps, 1.0::R)
		     else 1.0::R
	ncpL:(Integer, R)->R :=
	     x:I := #1
	     alpha:R := #2
	     if x=lo then 0.0::R
	     else 
		p:R := pnhyper(x, 1, true)
		if p>alpha then
			f:(R->R) := pnhyper(x,#1,true) - alpha
			ridder(f, 0,1)
		else if p<alpha then
			f:(R->R) := pnhyper(x,1/#1,true) - alpha
			1/ridder(f, doubleEps,1.0::R)
		else 1.0::R
 	pValue:R :=
 	     if alternative="less" then pnhyper(x00, OR,false) 
	     else if alternative="greater" then pnhyper(x00, OR,true)
	     else if alternative="two-sided" then
		relErr:= 1+1.0e-7::R
		dn:= dnhyper(OR)
		dstar:= dn(x00-lo+1)*relErr
		reduce("+",[di for di in dn | di<dstar])
	     else -1.0::R
	cInterval:List R :=
	       if confInt then 
		   if alternative="less" then  [0.0::R, ncpU(x00, 1.0::R-confLevel)]
		   else if alternative="greater" then [ncpL(x00, 1.0::R-confLevel), plusInfinity]
		   else if alternative="two-sided" then 
		      alpha:=(1-confLevel)/2
		      [ncpL(x00, alpha), ncpU(x00, alpha)]
		else [-1.0::R,-1.0::R]
		--else []
  	estimate:= mle(x00)
	[pValue, cInterval, estimate]
   testTolerance(x, y, atol) ==
    if abs(x-y) <= atol then true else false
   test1() ==
    testTolerance(2*ridder(cos,0.0::R,2.0::R),pi()$Pi::R, 1.0e-18)
   test2() == testTolerance(choose(100, 5)::R, 75287520::R, 0)
   test3() == testTolerance(dhyper(5, 10, 7, 8)::R, 0.3628137::R, 1.0e-7)
   test4() == testTolerance(log(dhyper(5, 10, 7, 8)::R),-1.013866::R,  1.0e-7)
   test5() == testTolerance(phyper(5, 10, 7, 8, true)::R,0.7821884::R, 1.0e-7)
   test6() == testTolerance(phyper(5, 10, 7, 8, false)::R,0.2178116::R, 1.0e-7) 
   test7() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).PValue,
       0.2575, 1.0e-3)
   test8() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).CI.1,
       0.5383996, 1.0e-6)
   test9() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).CI.2,
       7.4363242, 1.0e-4)
   test10() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).Estimate,
       1.971640, 1.0e-4)
   alltests() == [test1(), test2(), test3(), test4(), test5(), test6(),
       test7(), test8(), test9(), test10()]
\end{spad}

Using this code in Axiom:

\begin{axiom}
-- test code is correct
alltests()
-- show the example
fisherTest(10,10,10,20,"two-sided",1.0,true,0.95)
\end{axiom}

The Boot translation was more fiddly - but, then again, I had never used Boot before.

\begin{boot}
doubleFloat(x) == COERCE(x,'DOUBLE_-FLOAT)
DF(x) == COERCE(x,'DOUBLE_-FLOAT)
ridder(func, x1, x2) ==
    --x2:=DF(x2)
    eps:= DF(1.0e-16)
    maxit:= 30
    fl := DF(FUNCALL(func,x1))
    fh := DF(FUNCALL(func,x2))
    xl := x1
    xh := x2
    ans := DF(-1.11e20)
    xnew := 0.0e0
    iterNum:= 0
    if fl=0.0 then return x1
    else if fh=0.0 then return x2
    else if (fl*fh) > 0.0 then error "Initial points are not either side of zero."
    --if (fl*fh) < 0.0 then
    else repeat
		xm:= 0.5 *(xl+xh)
		fm:= FUNCALL(func,xm)
		ss:= SQRT((fm*fm) - (fl*fh))
		if ss =0.0 then return ans
		xnew:= xm + (((xm - xl) * (if (fl>fh) then 1.0 else -1.0) * fm) / ss)
		if ABS(xnew-ans) <= eps then return ans
		ans:= xnew 
		fnew:= DF(FUNCALL(func,ans))
		if fnew=0.0 then return ans
		if msign(fm,fnew) ^= fm then
		    xl:= xm 
		    fl:= fm 
		    xh:= ans 
		    fh:= fnew
		else if msign(fl, fnew) ^= fl then
		    xh:= ans 
		    fh:= fnew
		else if msign(fh, fnew) ^= fh then
		    xl:= ans 
		    fl:= fnew
		iterNum:=iterNum+1
		if iterNum >=maxit then 
			error "Maximum iterations exceeded"
		--if verbose then FORMAT(true,"~,8f ~,8f ~,8f ~,8f~%", xl, xh, fl, fh)$Lisp
		if ABS(xh-xl) <= eps then return ans
msign(x, y) ==
    (ABS x) * (if y>0.0 then 1.0 else if y<0.0 then -1.0 else 0.0)
choose(n, x) ==
    total := 1
    for denom in 1..x repeat
	total:=total*(n-denom+1)/denom
    return total
   --chooseNew(n, x) == product((n-i+1)::Fraction Integer/i::Fraction Integer,i=1..x)
dhyper(x, m, n, k) ==
    DF(choose(m, x) * choose(n, k - x)) / choose(m + n, k)
-- reduce(func,list) == 
-- 	value := list.0
-- 	for i in 1..(#list-1) repeat
-- 		value:=FUNCALL(func,value,list.i))
-- 	value
phyper(x, m, n, k, lowerTail) ==
    --total:Fraction Integer:=0/1
    if lowerTail then 
	+/[dhyper(i, m, n, k) for i in 1..x]
    else 
	+/[dhyper(i, m, n, k) for i in (x+1)..k]
dnhyper(ncp,logdc,support) ==
     d := [DF(logdc.i+LOG(ncp)*support.i) for i in 0..(#logdc-1)]
     maxd := APPLY(FUNCTION(MAX),d)
     d2 :=[EXP(di-maxd) for di in d]
     sumd2 := +/d2
     [d2i/sumd2 for d2i in d2]
testTolerance(x, y, atol) ==
    if ABS(x-y) <= atol then true else false
test1() == testTolerance(2*ridder('COS,0.0,2.0),3.1415926535897932385, 1.0e-7)
test2() == testTolerance(choose(100, 5), 75287520, 0)
test3() == testTolerance(dhyper(5, 10, 7, 8), 0.3628137, 1.0e-7)
test4() == testTolerance(LOG(dhyper(5, 10, 7, 8)),-1.013866,  1.0e-7)
test5() == testTolerance(phyper(5, 10, 7, 8, true),0.7821884, 1.0e-7)
test6() == testTolerance(phyper(5, 10, 7, 8, false),0.2178116, 1.0e-7) 
fisherTest(a,b,c,d, alternative, OR, confInt, confLevel) == main where
  main() ==
	$m := a+c -- first column
	$n := b+d -- second column
	$k := a+b -- first row
	$x00 := a
	$lo := MAX(0, $k-$n)
	$hi := MIN($k, $m)
	$support := [i for i in $lo..$hi]
	$logdc := [LOG(dhyper(i, $m, $n, $k)) for i in $support]
	$doubleEps := 1.0e-10
	$plusInfinity := 1.0e10
	pvalue :=
		if alternative='"less" then pnhyper($x00, OR,false) 
		else if alternative='"greater" then pnhyper($x00, OR,true)
		else if alternative='"two-sided" then
			relErr:= 1+1.0e-7
			d:= dnhyper(OR,$logdc,$support)
			dstar:= ELT(d,$x00-$lo)*relErr
			+/[di for di in d | di<dstar]
		else -1.0 -- no match
	estimate :=
		if $x00=$lo then 0
		-- else if $x00=hi then return($plusInfinity)
		else
			mu:= mnhyper(1)
			if mu>$x00 then ridder(FUNCTION(f1),0,1) 
			else if mu<$x00 then 1/ridder(FUNCTION(f2),$doubleEps,1) 
			else 1
 	interval :=
 	       if confInt then 
 			$alpha := 1 - confLevel
			if alternative='"less" then [0, ncpU($x00)]
			else if alternative='"greater" then [ncpL($x00), $plusInfinity]
			else if alternative='"two-sided" then 
				$alpha :=(1-confLevel)/2.0
				[ncpL($x00), ncpU($x00)]
			else [-1,-1] 
 		else [-2,-2]
	[pvalue,interval,estimate]
  pnhyper (q,ncp,upperTail) ==
	     if ncp=1 then 
		if upperTail then phyper(q-1, $m, $n, $k, false)
		else phyper(q, $m, $n, $k, true)
	     else if ncp=0 then
		if upperTail then 
		    if q<=$lo then 1 else 0
		else if q>=$lo then 1 else 0
-- 	     else if ncp=$plusInfinity then
-- 		if upperTail then 
-- 		    if q<=hi then return(1) else return(0)
-- 		else if q>= hi then return(1) else return(0)
	     else 
		d:= dnhyper(ncp, $logdc, $support)
		if upperTail then
		      +/[d.i for i in 0..(#d-1) | $support.i>=q]
		else +/[d.i for i in 0..(#d-1) | $support.i<=q]
  mnhyper(ncp) ==
	     if ncp=0.0 then $lo
 	     --if ncp=$plusInfinity then return(hi::R)
	     else
		d := dnhyper(ncp,$logdc,$support)
		+/[si*di for di in d for si in $support]
  f1(u) == mnhyper(u) - $x00
  f2(u) == mnhyper(1/u) - $x00
  ncpU x ==
	     --if x=$hi then $plusInfinity
	     p:= pnhyper(x, 1.0, false)
	     if p<$alpha then 
		ridder(FUNCTION(fu1),0.0,1.0)
	     else if p>$alpha then 
		1/ridder(FUNCTION(fu2), $doubleEps,1)
	     else 1
  fu1 u == pnhyper($x00,u,false) - $alpha
  fu2 u == pnhyper($x00,1/u,false) - $alpha
  ncpL x == 
	     if x=$lo then 0
	     else 
		p:= pnhyper(x, 1, true)
		if p>$alpha then ridder(FUNCTION(fl1), 0,1)
		else if p<$alpha then 1/ridder(FUNCTION(fl2), $doubleEps,1)
		else 1
  fl1 u == pnhyper($x00,u,true) - $alpha
  fl2 u == pnhyper($x00,1/u,true) - $alpha
test7() == fisherTest(10,10,10,20,'"two-sided",1,true,0.95)
alltests() == [test1(), test2(), test3(), test4(), test5(), test6()]
\end{boot}

Using this code in Axiom:

\begin{axiom}
alltests()$Lisp
fisherTest(10,10,10,20,"two-sided",1,true,0.95::SF)$Lisp
\end{axiom}

The implementation in Common Lisp was a more direct translation of the R code:

\begin{lisp}
;; from cl-statistics.lisp
(defun safe-exp (x)
  "Eliminates floating point underflow for the exponential function.
Instead, it just returns 0.0d0"
  (setf x (coerce x 'double-float))
  (if (< x (log least-positive-double-float))
      0.0d0
      (exp x)))
(defun ridder (func x1 x2 &key (eps 1.0d-16) (maxit 30) (verbose nil))
  (let (
	(fl (funcall func x1))
	(fh (funcall func x2))
	(xl x1)
	(xh x2)
	(ans -1.11d30)
	(xnew 0.0d0)
	(iter-num 0)
	)
    (cond
     ((= fl 0) x1)
     ((= fh 0) x2)
     ((> (* fl fh) 0.0d0) 
      (error "Functions of the start points are not either side of zero."))
     ((< (* fl fh) 0.0d0) 
      (loop
       (let* (
	      (xm (* 0.5d0 (+ xl xh)))
	      (fm (funcall func xm))
	      (ss (sqrt (- (* fm fm) (* fl fh))))
	      )
	 (if (= ss 0.0d0) (return ans))
	 (setf xnew (+ xm (/ (* (- xm xl) (if (> fl fh) 1.0d0 -1.0d0) fm) ss)))
	 (if (<= (abs (- xnew ans)) eps) (return ans))
	 (setf ans xnew fnew (funcall func ans))
	 (if (= fnew 0.0d0) (return ans))
	 (cond ((not (= (msign fm fnew) fm))
		(setf xl xm fl fm xh ans fh fnew))
	       ((not (= (msign fl fnew) fl))
		(setf xh ans fh fnew))
	       ((not (= (msign fh fnew) fh))
		(setf xl ans fl fnew)))
	 (incf iter-num)
	 (if (>= iter-num maxit) 
	     (return (values nil "Maximum iterations exceeded"))) ;; (error)?
	 (if verbose (format t "~,8f ~,8f ~,8f ~,8f~%" xl xh fl fh))
	 (if (<= (abs (- xh xl)) eps) (return ans))))))))
(defun msign (x y)
  (* (abs x) (cond ((> y 0.0d0) 1.0d0) ((< y 0.0d0) -1.0d0) (t 0.0d0))))
;;(- (* (ridder #'cos 0.0d0 2.0d0) 2.0d0) pi)

(defun choose (n x)
  (loop 
     for denom from 1 to x 
     and numerator from n downto (- n (1- x)) 
     and total = 1 then (* total (/ numerator denom)) 
     finally (return total)))
(defun dhyper (x m n k &key (log nil))
  (let ((val
	 (/ (* (choose m x) (choose n (- k x))) (choose (+ m n) k))))
    (if log 
	(log (coerce val 'double-float)) 
	val)))
(defun phyper (x m n k &key (lower-tail t))
  (if lower-tail 
      (loop for i from 1 to x summing (dhyper i m n k))
      (loop for i from (1+ x) to k summing (dhyper i m n k))))

(defun fisher-test (x &key (alternative 'two-sided) (or 1.0d0)
		    (conf-int t)
		    (conf-level 0.95d0) (uniroot #'ridder))
    "Fisher's exact test for a 2x2 integer array.
This is a hand translation of R's fisher.test() 
making use of CL's large integers for the hypergeometric distribution"
    (let* ((m (loop for i upto 1 summing (aref x i 0)))
	   (n (loop for i upto 1 summing (aref x i 1)))
	   (k (loop for i upto 1 summing (aref x 0 i)))
	   (x00 (aref x 0 0)) ; cf replacing x by (aref x 0 0)
	   (lo (max 0 (- k n)))
	   (hi (min k m))
	   (support (loop for i from lo to hi collect i))
	   (log-dc (loop for i in support 
		      collect (dhyper i m n k :log t)))
	   (double-eps 1.0d-50))
      (labels 
	  ((dnhyper (ncp)
	     (setf ncp (coerce ncp 'double-float))
	     (let* ((d (loop for i in log-dc and j in support
			  collect (+ i (* (log  ncp) j))))
		    (max-d (apply #'max d))
		    (d2 (loop for i in d collect 
			     (safe-exp (- i max-d)))) ;; NB: safe-exp used here
		    (sum-d2 (reduce #'+ d2)))
	       (loop for i in d2 collect (/ i sum-d2))))
	   (mnhyper (ncp)
	     (cond
	       ((= ncp 0) lo)
	       ((equal ncp 'infinity) hi)
	       (t (loop for i in support and j in (dnhyper ncp)
		       summing (* i j)))))
	   (pnhyper (q ncp &key (upper-tail nil))
	     (cond 
	       ((= ncp 1)
		(if upper-tail 
		    (coerce (phyper (1- x00) m n k :lower-tail nil) 'double-float)
		    (coerce (phyper x00 m n k) 'double-float)))
	       ((= ncp 0)
		(if upper-tail 
		    (if (<= q lo) 1 0)
		    (if (>= q lo) 1 0)))
	       ((equal ncp 'infinity)
		(if upper-tail 
		    (if (<= q hi) 1 0)
		    (if (>= q hi) 1 0)))
	       (t 
		(let ((d (dnhyper ncp)))
		  (if upper-tail
		      (loop for d-i in d and support-i in support
			   when (>= support-i q)
			   summing d-i)
		      (loop for d-i in d and support-i in support
			   when (<= support-i q)
			   summing d-i))))))
	   (mle (x)
	     (cond
	       ((= x lo) 0)
	       ((= x hi) 'infinity)
	       (t
		(let ((mu (mnhyper 1)))
		  (cond 
		    ((> mu x)
		     (funcall uniroot (lambda (u) (- (mnhyper u) x))
			      0 1))
		    ((< mu x)
		     (/ (funcall uniroot (lambda (u) (- (mnhyper (/ u)) x))
			      double-eps 1)))
		    (t 1))))))
	   (ncp-u (x alpha)
	     (and (= x hi) 'infinity)
	     (let ((p (pnhyper x 1)))
	       (cond 
		 ((< p alpha)
		  (funcall uniroot (lambda (u) (- (pnhyper x u) alpha))
			 0 1))
		 ((> p alpha)
		  (/ (funcall uniroot (lambda (u) (- (pnhyper x (/ u)) alpha))
			    double-eps 1)))
		 (t 1))))
	   (ncp-l (x alpha)
	     (and (= x lo) 0)
	     (let ((p (pnhyper x 1 :upper-tail t)))
	       (cond 
		 ((> p alpha)
		  (funcall uniroot (lambda (u) (- (pnhyper x u :upper-tail t) alpha))
			 0 1))
		 ((< p alpha)
		  (/ (funcall uniroot (lambda (u) (- (pnhyper x (/ u) :upper-tail t) alpha))
			    double-eps 1)))
		 (t 1)))))
	(let ((p-value 
	       (ecase alternative
		 (less (pnhyper x00 or))
		 (greater (pnhyper x00 or :upper-tail t)) 
		 (two-sided 
		  (let* ((relErr  (1+ 1.0d-7))
			 (d (dnhyper or))
			 (dstar (* (elt d (- x00 lo)) relErr)))
		    (loop for di in d when (< di dstar) summing di)))))
	      (c-interval
	       (if conf-int
		   (ecase alternative
		     (less (list 0 (ncp-u x00 (- 1 conf-level))))
		     (greater (list (ncp-l x00 (- 1 conf-level)) 'infinity))
		     (two-sided 
		      (let ((alpha (/ (- 1 conf-level) 2)))
			(list (ncp-l x00 alpha) (ncp-u x00 alpha)))))
		   nil))
	      (estimate (mle x00)))
	  (values p-value c-interval estimate)))))
;;(fisher-test #2a((10 10) (10 20)))
\end{lisp}

With output:

\begin{verbatim}
CL-USER> (fisher-test #2a((10 10) (10 20)))
0.2575492428109829d0
(0.5383993816781727d0 7.4363408387439875d0)
1.9716269432603386d0
\end{verbatim}

How useful are the different CAS languages for implementing numerical routines? Prompted by a comparison of R and C for implementing Fisher's exact test for 2x2 tables (http://fluff.info/blog/arch/00000172.htm), I thought that it would be interesting to implement this particular test in Spad, Boot and Common Lisp (see below). Each set of code was required to implement a univariate root finder and the hypergeometric distribution to calculate the p-value under different alternatives, together with the 95% confidence interval and the maximum likelihood estimator for the odds ratio. The reference implementation is R, where the code and output would be:

LatexWiki Image

To summarise, all three languages (Spad, Boot and Common Lisp) provide arbitrary length integers and fractions, ensuring that the hypergeometric distribution was straightforward to implement. The Lisp and R implementations were very similar, which is not surprising, given that the two languages are closely related. In contrast, the nested functions in Spad seemed clunky, with the requirement to use the #1 and #2 argument references, although I did appreciate Spad's lexical scoping and facility to fall back to an algebraic analysis. An Aldor implementation would be cleaner than the Spad implementation. Boot's implementation was initially difficult, as I was unclear how to pass values to the nested functions with using function argument (which was not possible for the univariate root finder :-(). The use of the "$" prefix on derived variables seemed clunky - and I was surprised to find that function arguments were lexically scoped. Finally, my version of Fricas:Boot defaulted to single precision floats, which caused problems with precision. Type specification in OpenAxiom:Boot and more recent versions of Fricas would negate the need for explicit coercion to double-floats.

This begs the question: when would one use any of these languages for mixed numerical/algebraic analysis? First, Boot is the least likely to be used, although it does play closely and well with Common Lisp (a la Maxima). One could code for numerical analysis in Boot and Common Lisp - however Boot's lack of lexical scoping may be a detraction. Moreover, by my understanding, Lisp or Boot are unable to evaluate Spad or Axiom functions. Second, for an R user, Spad's type system seems fussy and debugging can be slow. The lack of ability to call Spad functions from Lisp or Boot is a major restriction, requiring that a large body of code in Lisp (or Fortran via f2cl), such as for optimisation, would need to be hand translated to Spad (see [SandBoxMLE]? for an example). As a final note, Maxima would seem to address many of these concerns - although it lacks an elegant type system and I haven't learnt to appreciate its scoping rules :-).

spad

)abbrev package TESTP TestPackage
R ==> Float
I ==> Integer
fisherRec ==> Record(PValue:R, CI:List R, Estimate:R)
TestPackage: with
   ridder: (R->R,R,R) -> R
   msign: (R,R) -> R
   choose:(I, I)->Fraction I 
   --chooseNew:(Integer, Integer)->Fraction Integer 
   dhyper:(I, I, I, I)->Fraction Integer
   phyper:(I, I, I, I, Boolean)->Fraction Integer
   fisherTest:(I,I,I,I, String, R, Boolean, R)-> fisherRec
   testTolerance:(R, R, R)->Boolean   
   test1: () -> Boolean
   test2: () -> Boolean
   test3: () -> Boolean
   test4: () -> Boolean
   test5: () -> Boolean
   test6: () -> Boolean
   test7: () -> Boolean
   test8: () -> Boolean
   test9: () -> Boolean
   test10: () -> Boolean
   alltests: () -> List Boolean
  == add
   import TrigonometricFunctionCategory -- for test1()
   --import OrderedCompletion(Float) for plusInfinity()$OrderedCompletion(Float)
   ridder(func, x1, x2) ==
    eps:= 1.0e-16::R
    maxit:= 30::Integer
    --verbose:= false
    fl:R := func x1
    fh:R := func x2
    xl:R := x1
    xh:R := x2
    ans:R := -1.11e30::R
    xnew:R := 0.0e0::R
    iterNum:= 0::Integer
    if fl=0.0::R then return x1
    else if fh=0.0::R then return x2
    else if (fl*fh) > 0.0::R then error "Initial points are not either side of zero."
    --if (fl*fh) < 0.0 then
    else repeat
                xm:= 0.5::R *(xl+xh)
                fm:= func xm
                ss:= sqrt((fm*fm) - (fl*fh))
                if ss =0.0::R then return ans
                xnew:= xm + (((xm - xl) * (if (fl>fh) then 1.0::R else -1.0::R) * fm) / ss)
                if abs(xnew-ans) <= eps then return ans
                ans:= xnew 
                fnew:= func ans
                if fnew=0.0::R then return ans
                if msign(fm,fnew) ~= fm then
                    xl:= xm 
                    fl:= fm 
                    xh:= ans 
                    fh:= fnew
                else if msign(fl, fnew) ~= fl then
                    xh:= ans 
                    fh:= fnew
                else if msign(fh, fnew) ~= fh then
                    xl:= ans 
                    fl:= fnew
                iterNum:=iterNum+1::Integer
                if iterNum >=maxit then 
                        error "Maximum iterations exceeded"
                --if verbose then FORMAT(true,"~,8f ~,8f ~,8f ~,8f~%", xl, xh, fl, fh)$Lisp
                if abs(xh-xl) <= eps then return ans
   msign(x, y) ==
    (abs x) * (if y>0.0::R then 1.0::R else if y<0.0::R then -1.0::R else 0.0::R)
   choose(n, x) ==
    total:Fraction Integer := 1/1
    for denom in 1..x repeat
        total:=total*((n-denom+1)/(denom))::Fraction Integer
    return total
   --chooseNew(n, x) == product((n-i+1)::Fraction Integer/i::Fraction Integer,i=1..x)
   dhyper(x, m, n, k) ==
    choose(m, x) * choose(n, k - x) / choose(m + n, k)
   phyper(x, m, n, k, lowerTail) ==
    i:PositiveInteger
    --total:Fraction Integer:=0/1
    if lowerTail then 
        reduce("+",[dhyper(i, m, n, k) for i in 1..x])
    else 
        reduce("+",[dhyper(i, m, n, k) for i in (x+1)..k])
   fisherTest(a,b,c,d, alternative, OR, confInt, confLevel) ==
        m:I := a+c -- first column
        n:I := b+d -- second column
        k:I := a+b -- first row
        x00:I := a
        lo:I := max(0, k-n)
        hi:I := min(k, m)
        support:List I := [i for i in lo..hi]
        logdc:List R:= [log(dhyper(i, m, n, k)::R) for i in support]
        doubleEps:R := 1.0e-50::R 
        plusInfinity:R := 1.0e6::R -- arbitrary
        dnhyper:(R->List R) :=
             ncp:R := #1
             d:List R := [logdc(i)+log(ncp)*support(i)::R for i in 1..#logdc]
             maxd:R := reduce(max,d)
             d2:List R :=[exp(di-maxd) for di in d]
             sumd2:R := reduce("+",d2)
             [d2i/sumd2 for d2i in d2]
        mnhyper:(R->R) :=
             ncp:R := #1
             if ncp=0.0::R then lo::R
             --else if ncp=%plusInfinity then hi::R
             else 
                d:List R := dnhyper(ncp)
                reduce("+",[support(i)::R*d(i) for i in 1..#d])
        pnhyper:((Integer,R,Boolean)->R) :=
             q:I := #1
             ncp:R := #2
             upperTail:Boolean := #3
             if ncp=1.0 then 
                if upperTail then phyper(q-1, m, n, k, false)::R
                else phyper(q, m, n, k, true)::R
             else if ncp=0.0 then
                if upperTail then 
                    if q<=lo then 1.0::R else 0.0::R
                else if q>=lo then 1.0::R else 0.0::R
--           else if ncp=%plusInfinity then
--              if upperTail then 
--                  if q<=hi then 1.0::R else 0.0::R
--              else if q>= hi then 1.0::R else 0.0::R
             else 
                d:List R := dnhyper(ncp)
                if upperTail then
                      reduce("+",[d(i) for i in 1..#d | support(i)>=q])
                else reduce("+",[d(i) for i in 1..#d | support(i)<=q])
        mle:(I->R) :=
             x:I := #1
             if x=lo then 0.0::R
             else if x=hi then plusInfinity
             else
                mu:R := mnhyper(1.0::R)
                if mu>x::R then 
                        f:(R->R) := mnhyper(#1) - x::R
                        ridder(f,0,1)
                else if mu<x::R then
                        f:(R->R) := mnhyper(1/#1) - x::R
                        1/ridder(f,doubleEps,1.0::R)
                else 1.0::R
        ncpU:(I,R)->R :=
             x:I := #1
             alpha:R := #2
             if x=hi then plusInfinity
             else
                     p:R := pnhyper(x, 1.0::R, false)
                     if p<alpha then 
                        f:(R->R) := pnhyper(x,#1,false) - alpha
                        ridder(f, 0.0::R, 1.0::R)
                     else if p>alpha then 
                                f:(R->R) := pnhyper(x,1/#1,false) - alpha
                                1/ridder(f, doubleEps, 1.0::R)
                     else 1.0::R
        ncpL:(Integer, R)->R :=
             x:I := #1
             alpha:R := #2
             if x=lo then 0.0::R
             else 
                p:R := pnhyper(x, 1, true)
                if p>alpha then
                        f:(R->R) := pnhyper(x,#1,true) - alpha
                        ridder(f, 0,1)
                else if p<alpha then
                        f:(R->R) := pnhyper(x,1/#1,true) - alpha
                        1/ridder(f, doubleEps,1.0::R)
                else 1.0::R
        pValue:R :=
             if alternative="less" then pnhyper(x00, OR,false) 
             else if alternative="greater" then pnhyper(x00, OR,true)
             else if alternative="two-sided" then
                relErr:= 1+1.0e-7::R
                dn:= dnhyper(OR)
                dstar:= dn(x00-lo+1)*relErr
                reduce("+",[di for di in dn | di<dstar])
             else -1.0::R
        cInterval:List R :=
               if confInt then 
                   if alternative="less" then  [0.0::R, ncpU(x00, 1.0::R-confLevel)]
                   else if alternative="greater" then [ncpL(x00, 1.0::R-confLevel), plusInfinity]
                   else if alternative="two-sided" then 
                      alpha:=(1-confLevel)/2
                      [ncpL(x00, alpha), ncpU(x00, alpha)]
                else [-1.0::R,-1.0::R]
                --else []
        estimate:= mle(x00)
        [pValue, cInterval, estimate]
   testTolerance(x, y, atol) ==
    if abs(x-y) <= atol then true else false
   test1() ==
    testTolerance(2*ridder(cos,0.0::R,2.0::R),pi()$Pi::R, 1.0e-18)
   test2() == testTolerance(choose(100, 5)::R, 75287520::R, 0)
   test3() == testTolerance(dhyper(5, 10, 7, 8)::R, 0.3628137::R, 1.0e-7)
   test4() == testTolerance(log(dhyper(5, 10, 7, 8)::R),-1.013866::R,  1.0e-7)
   test5() == testTolerance(phyper(5, 10, 7, 8, true)::R,0.7821884::R, 1.0e-7)
   test6() == testTolerance(phyper(5, 10, 7, 8, false)::R,0.2178116::R, 1.0e-7) 
   test7() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).PValue,
       0.2575, 1.0e-3)
   test8() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).CI.1,
       0.5383996, 1.0e-6)
   test9() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).CI.2,
       7.4363242, 1.0e-4)
   test10() == testTolerance(fisherTest(10,10,10,20,"two-sided",1.0,true,0.95).Estimate,
       1.971640, 1.0e-4)
   alltests() == [test1(), test2(), test3(), test4(), test5(), test6(),
       test7(), test8(), test9(), test10()]

spad

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/303933172962875058-25px001.spad using 
      old system compiler.
   TESTP abbreviates package TestPackage 
   processing macro definition R ==> Float 

   processing macro definition I ==> Integer 

   processing macro definition fisherRec ==> Record(PValue: R,CI: List R,Estimate: R) 

------------------------------------------------------------------------
   initializing NRLIB TESTP for TestPackage 
   compiling into NRLIB TESTP 
   importing TrigonometricFunctionCategory
   compiling exported ridder : (Float -> Float,Float,Float) -> Float
Time: 0.06 SEC.

   compiling exported msign : (Float,Float) -> Float
Time: 0.02 SEC.

   compiling exported choose : (Integer,Integer) -> Fraction Integer
Time: 0.04 SEC.

   compiling exported dhyper : (Integer,Integer,Integer,Integer) -> Fraction Integer
Time: 0.01 SEC.

   compiling exported phyper : (Integer,Integer,Integer,Integer,Boolean) -> Fraction Integer
Time: 0.03 SEC.

   compiling exported fisherTest : (Integer,Integer,Integer,Integer,String,Float,Boolean,Float) -> Record(PValue: Float,CI: List Float,Estimate: Float)
Time: 0.20 SEC.

   compiling exported testTolerance : (Float,Float,Float) -> Boolean
Time: 0.01 SEC.

   compiling exported test1 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test2 : () -> Boolean
Time: 0 SEC.

   compiling exported test3 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test4 : () -> Boolean
Time: 0.05 SEC.

   compiling exported test5 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test6 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test7 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test8 : () -> Boolean
Time: 0.01 SEC.

   compiling exported test9 : () -> Boolean
Time: 0.02 SEC.

   compiling exported test10 : () -> Boolean
Time: 0 SEC.

   compiling exported alltests : () -> List Boolean
Time: 0.01 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |TestPackage| REDEFINED

;;;     ***       |TestPackage| REDEFINED
Time: 0 SEC.

   Warnings: 
      [1] ridder:  xh has no value
      [2] ridder:  xl has no value
      [3] fisherTest: The conditional modes (List (Float)) and (Integer) conflict

   Cumulative Statistics for Constructor TestPackage
      Time: 0.51 seconds

   finalizing NRLIB TESTP 
   Processing TestPackage for Browser database:
--->-->TestPackage((ridder (R (Mapping R R) R R))): Not documented!!!!
--->-->TestPackage((msign (R R R))): Not documented!!!!
--->-->TestPackage((choose ((Fraction I) I I))): Not documented!!!!
--->-->TestPackage((dhyper ((Fraction (Integer)) I I I I))): Not documented!!!!
--->-->TestPackage((phyper ((Fraction (Integer)) I I I I (Boolean)))): Not documented!!!!
--->-->TestPackage((fisherTest (fisherRec I I I I (String) R (Boolean) R))): Not documented!!!!
--->-->TestPackage((testTolerance ((Boolean) R R R))): Not documented!!!!
--->-->TestPackage((test1 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test2 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test3 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test4 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test5 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test6 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test7 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test8 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test9 ((Boolean)))): Not documented!!!!
--->-->TestPackage((test10 ((Boolean)))): Not documented!!!!
--->-->TestPackage((alltests ((List (Boolean))))): Not documented!!!!
--->-->TestPackage(constructor): Not documented!!!!
--->-->TestPackage(): Missing Description
------------------------------------------------------------------------
   TestPackage is now explicitly exposed in frame initial 
   TestPackage will be automatically loaded when needed from 
      /var/zope2/var/LatexWiki/TESTP.NRLIB/code

Using this code in Axiom:

axiom

-- test code is correct
alltests()

(1)

Type: List Boolean

axiom

-- show the example
fisherTest(10,10,10,20,"two-sided",1.0,true,0.95)

(2)

Type: Record(PValue?: Float,CI: List Float,Estimate: Float)

The Boot translation was more fiddly - but, then again, I had never used Boot before.

boot

doubleFloat(x) == COERCE(x,'DOUBLE_-FLOAT)
DF(x) == COERCE(x,'DOUBLE_-FLOAT)
ridder(func, x1, x2) ==
    --x2:=DF(x2)
    eps:= DF(1.0e-16)
    maxit:= 30
    fl := DF(FUNCALL(func,x1))
    fh := DF(FUNCALL(func,x2))
    xl := x1
    xh := x2
    ans := DF(-1.11e20)
    xnew := 0.0e0
    iterNum:= 0
    if fl=0.0 then return x1
    else if fh=0.0 then return x2
    else if (fl*fh) > 0.0 then error "Initial points are not either side of zero."
    --if (fl*fh) < 0.0 then
    else repeat
                xm:= 0.5 *(xl+xh)
                fm:= FUNCALL(func,xm)
                ss:= SQRT((fm*fm) - (fl*fh))
                if ss =0.0 then return ans
                xnew:= xm + (((xm - xl) * (if (fl>fh) then 1.0 else -1.0) * fm) / ss)
                if ABS(xnew-ans) <= eps then return ans
                ans:= xnew 
                fnew:= DF(FUNCALL(func,ans))
                if fnew=0.0 then return ans
                if msign(fm,fnew) ^= fm then
                    xl:= xm 
                    fl:= fm 
                    xh:= ans 
                    fh:= fnew
                else if msign(fl, fnew) ^= fl then
                    xh:= ans 
                    fh:= fnew
                else if msign(fh, fnew) ^= fh then
                    xl:= ans 
                    fl:= fnew
                iterNum:=iterNum+1
                if iterNum >=maxit then 
                        error "Maximum iterations exceeded"
                --if verbose then FORMAT(true,"~,8f ~,8f ~,8f ~,8f~%", xl, xh, fl, fh)$Lisp
                if ABS(xh-xl) <= eps then return ans
msign(x, y) ==
    (ABS x) * (if y>0.0 then 1.0 else if y<0.0 then -1.0 else 0.0)
choose(n, x) ==
    total := 1
    for denom in 1..x repeat
        total:=total*(n-denom+1)/denom
    return total
   --chooseNew(n, x) == product((n-i+1)::Fraction Integer/i::Fraction Integer,i=1..x)
dhyper(x, m, n, k) ==
    DF(choose(m, x) * choose(n, k - x)) / choose(m + n, k)
-- reduce(func,list) == 
--      value := list.0
--      for i in 1..(#list-1) repeat
--              value:=FUNCALL(func,value,list.i))
--      value
phyper(x, m, n, k, lowerTail) ==
    --total:Fraction Integer:=0/1
    if lowerTail then 
        +/[dhyper(i, m, n, k) for i in 1..x]
    else 
        +/[dhyper(i, m, n, k) for i in (x+1)..k]
dnhyper(ncp,logdc,support) ==
     d := [DF(logdc.i+LOG(ncp)*support.i) for i in 0..(#logdc-1)]
     maxd := APPLY(FUNCTION(MAX),d)
     d2 :=[EXP(di-maxd) for di in d]
     sumd2 := +/d2
     [d2i/sumd2 for d2i in d2]
testTolerance(x, y, atol) ==
    if ABS(x-y) <= atol then true else false
test1() == testTolerance(2*ridder('COS,0.0,2.0),3.1415926535897932385, 1.0e-7)
test2() == testTolerance(choose(100, 5), 75287520, 0)
test3() == testTolerance(dhyper(5, 10, 7, 8), 0.3628137, 1.0e-7)
test4() == testTolerance(LOG(dhyper(5, 10, 7, 8)),-1.013866,  1.0e-7)
test5() == testTolerance(phyper(5, 10, 7, 8, true),0.7821884, 1.0e-7)
test6() == testTolerance(phyper(5, 10, 7, 8, false),0.2178116, 1.0e-7) 
fisherTest(a,b,c,d, alternative, OR, confInt, confLevel) == main where
  main() ==
        $m := a+c -- first column
        $n := b+d -- second column
        $k := a+b -- first row
        $x00 := a
        $lo := MAX(0, $k-$n)
        $hi := MIN($k, $m)
        $support := [i for i in $lo..$hi]
        $logdc := [LOG(dhyper(i, $m, $n, $k)) for i in $support]
        $doubleEps := 1.0e-10
        $plusInfinity := 1.0e10
        pvalue :=
                if alternative='"less" then pnhyper($x00, OR,false) 
                else if alternative='"greater" then pnhyper($x00, OR,true)
                else if alternative='"two-sided" then
                        relErr:= 1+1.0e-7
                        d:= dnhyper(OR,$logdc,$support)
                        dstar:= ELT(d,$x00-$lo)*relErr
                        +/[di for di in d | di<dstar]
                else -1.0 -- no match
        estimate :=
                if $x00=$lo then 0
                -- else if $x00=hi then return($plusInfinity)
                else
                        mu:= mnhyper(1)
                        if mu>$x00 then ridder(FUNCTION(f1),0,1) 
                        else if mu<$x00 then 1/ridder(FUNCTION(f2),$doubleEps,1) 
                        else 1
        interval :=
               if confInt then 
                        $alpha := 1 - confLevel
                        if alternative='"less" then [0, ncpU($x00)]
                        else if alternative='"greater" then [ncpL($x00), $plusInfinity]
                        else if alternative='"two-sided" then 
                                $alpha :=(1-confLevel)/2.0
                                [ncpL($x00), ncpU($x00)]
                        else [-1,-1] 
                else [-2,-2]
        [pvalue,interval,estimate]
  pnhyper (q,ncp,upperTail) ==
             if ncp=1 then 
                if upperTail then phyper(q-1, $m, $n, $k, false)
                else phyper(q, $m, $n, $k, true)
             else if ncp=0 then
                if upperTail then 
                    if q<=$lo then 1 else 0
                else if q>=$lo then 1 else 0
--           else if ncp=$plusInfinity then
--              if upperTail then 
--                  if q<=hi then return(1) else return(0)
--              else if q>= hi then return(1) else return(0)
             else 
                d:= dnhyper(ncp, $logdc, $support)
                if upperTail then
                      +/[d.i for i in 0..(#d-1) | $support.i>=q]
                else +/[d.i for i in 0..(#d-1) | $support.i<=q]
  mnhyper(ncp) ==
             if ncp=0.0 then $lo
             --if ncp=$plusInfinity then return(hi::R)
             else
                d := dnhyper(ncp,$logdc,$support)
                +/[si*di for di in d for si in $support]
  f1(u) == mnhyper(u) - $x00
  f2(u) == mnhyper(1/u) - $x00
  ncpU x ==
             --if x=$hi then $plusInfinity
             p:= pnhyper(x, 1.0, false)
             if p<$alpha then 
                ridder(FUNCTION(fu1),0.0,1.0)
             else if p>$alpha then 
                1/ridder(FUNCTION(fu2), $doubleEps,1)
             else 1
  fu1 u == pnhyper($x00,u,false) - $alpha
  fu2 u == pnhyper($x00,1/u,false) - $alpha
  ncpL x == 
             if x=$lo then 0
             else 
                p:= pnhyper(x, 1, true)
                if p>$alpha then ridder(FUNCTION(fl1), 0,1)
                else if p<$alpha then 1/ridder(FUNCTION(fl2), $doubleEps,1)
                else 1
  fl1 u == pnhyper($x00,u,true) - $alpha
  fl2 u == pnhyper($x00,1/u,true) - $alpha
test7() == fisherTest(10,10,10,20,'"two-sided",1,true,0.95)
alltests() == [test1(), test2(), test3(), test4(), test5(), test6()]

boot

Value = T
; (DEFUN |fisherTest,fl2| ...) is being compiled.
;; The variable |$x00| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$alpha| is undefined.
;; The compiler will assume this variable is a global.
; (DEFUN |fisherTest,ncpL| ...) is being compiled.
;; The variable |$lo| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$doubleEps| is undefined.
;; The compiler will assume this variable is a global.
; (DEFUN |fisherTest,mnhyper| ...) is being compiled.
;; The variable |$logdc| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$support| is undefined.
;; The compiler will assume this variable is a global.
; (DEFUN |fisherTest,pnhyper| ...) is being compiled.
;; The variable |$m| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$n| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$k| is undefined.
;; The compiler will assume this variable is a global.
; (DEFUN |fisherTest| ...) is being compiled.
;; The variable |$hi| is undefined.
;; The compiler will assume this variable is a global.
;; The variable |$plusInfinity| is undefined.
;; The compiler will assume this variable is a global.
Value = 17616

Using this code in Axiom:

axiom

alltests()$Lisp

(3)

Type: SExpression?

axiom

fisherTest(10,10,10,20,"two-sided",1,true,0.95::SF)$Lisp

(4)

Type: SExpression?

The implementation in Common Lisp was a more direct translation of the R code:

lisp

;; from cl-statistics.lisp
(defun safe-exp (x)
  "Eliminates floating point underflow for the exponential function.
Instead, it just returns 0.0d0"
  (setf x (coerce x 'double-float))
  (if (< x (log least-positive-double-float))
      0.0d0
      (exp x)))
(defun ridder (func x1 x2 &key (eps 1.0d-16) (maxit 30) (verbose nil))
  (let (
        (fl (funcall func x1))
        (fh (funcall func x2))
        (xl x1)
        (xh x2)
        (ans -1.11d30)
        (xnew 0.0d0)
        (iter-num 0)
        )
    (cond
     ((= fl 0) x1)
     ((= fh 0) x2)
     ((> (* fl fh) 0.0d0) 
      (error "Functions of the start points are not either side of zero."))
     ((< (* fl fh) 0.0d0) 
      (loop
       (let* (
              (xm (* 0.5d0 (+ xl xh)))
              (fm (funcall func xm))
              (ss (sqrt (- (* fm fm) (* fl fh))))
              )
         (if (= ss 0.0d0) (return ans))
         (setf xnew (+ xm (/ (* (- xm xl) (if (> fl fh) 1.0d0 -1.0d0) fm) ss)))
         (if (<= (abs (- xnew ans)) eps) (return ans))
         (setf ans xnew fnew (funcall func ans))
         (if (= fnew 0.0d0) (return ans))
         (cond ((not (= (msign fm fnew) fm))
                (setf xl xm fl fm xh ans fh fnew))
               ((not (= (msign fl fnew) fl))
                (setf xh ans fh fnew))
               ((not (= (msign fh fnew) fh))
                (setf xl ans fl fnew)))
         (incf iter-num)
         (if (>= iter-num maxit) 
             (return (values nil "Maximum iterations exceeded"))) ;; (error)?
         (if verbose (format t "~,8f ~,8f ~,8f ~,8f~%" xl xh fl fh))
         (if (<= (abs (- xh xl)) eps) (return ans))))))))
(defun msign (x y)
  (* (abs x) (cond ((> y 0.0d0) 1.0d0) ((< y 0.0d0) -1.0d0) (t 0.0d0))))
;;(- (* (ridder #'cos 0.0d0 2.0d0) 2.0d0) pi)

(defun choose (n x)
  (loop 
     for denom from 1 to x 
     and numerator from n downto (- n (1- x)) 
     and total = 1 then (* total (/ numerator denom)) 
     finally (return total)))
(defun dhyper (x m n k &key (log nil))
  (let ((val
         (/ (* (choose m x) (choose n (- k x))) (choose (+ m n) k))))
    (if log 
        (log (coerce val 'double-float)) 
        val)))
(defun phyper (x m n k &key (lower-tail t))
  (if lower-tail 
      (loop for i from 1 to x summing (dhyper i m n k))
      (loop for i from (1+ x) to k summing (dhyper i m n k))))

(defun fisher-test (x &key (alternative 'two-sided) (or 1.0d0)
                    (conf-int t)
                    (conf-level 0.95d0) (uniroot #'ridder))
    "Fisher's exact test for a 2x2 integer array.
This is a hand translation of R's fisher.test() 
making use of CL's large integers for the hypergeometric distribution"
    (let* ((m (loop for i upto 1 summing (aref x i 0)))
           (n (loop for i upto 1 summing (aref x i 1)))
           (k (loop for i upto 1 summing (aref x 0 i)))
           (x00 (aref x 0 0)) ; cf replacing x by (aref x 0 0)
           (lo (max 0 (- k n)))
           (hi (min k m))
           (support (loop for i from lo to hi collect i))
           (log-dc (loop for i in support 
                      collect (dhyper i m n k :log t)))
           (double-eps 1.0d-50))
      (labels 
          ((dnhyper (ncp)
             (setf ncp (coerce ncp 'double-float))
             (let* ((d (loop for i in log-dc and j in support
                          collect (+ i (* (log  ncp) j))))
                    (max-d (apply #'max d))
                    (d2 (loop for i in d collect 
                             (safe-exp (- i max-d)))) ;; NB: safe-exp used here
                    (sum-d2 (reduce #'+ d2)))
               (loop for i in d2 collect (/ i sum-d2))))
           (mnhyper (ncp)
             (cond
               ((= ncp 0) lo)
               ((equal ncp 'infinity) hi)
               (t (loop for i in support and j in (dnhyper ncp)
                       summing (* i j)))))
           (pnhyper (q ncp &key (upper-tail nil))
             (cond 
               ((= ncp 1)
                (if upper-tail 
                    (coerce (phyper (1- x00) m n k :lower-tail nil) 'double-float)
                    (coerce (phyper x00 m n k) 'double-float)))
               ((= ncp 0)
                (if upper-tail 
                    (if (<= q lo) 1 0)
                    (if (>= q lo) 1 0)))
               ((equal ncp 'infinity)
                (if upper-tail 
                    (if (<= q hi) 1 0)
                    (if (>= q hi) 1 0)))
               (t 
                (let ((d (dnhyper ncp)))
                  (if upper-tail
                      (loop for d-i in d and support-i in support
                           when (>= support-i q)
                           summing d-i)
                      (loop for d-i in d and support-i in support
                           when (<= support-i q)
                           summing d-i))))))
           (mle (x)
             (cond
               ((= x lo) 0)
               ((= x hi) 'infinity)
               (t
                (let ((mu (mnhyper 1)))
                  (cond 
                    ((> mu x)
                     (funcall uniroot (lambda (u) (- (mnhyper u) x))
                              0 1))
                    ((< mu x)
                     (/ (funcall uniroot (lambda (u) (- (mnhyper (/ u)) x))
                              double-eps 1)))
                    (t 1))))))
           (ncp-u (x alpha)
             (and (= x hi) 'infinity)
             (let ((p (pnhyper x 1)))
               (cond 
                 ((< p alpha)
                  (funcall uniroot (lambda (u) (- (pnhyper x u) alpha))
                         0 1))
                 ((> p alpha)
                  (/ (funcall uniroot (lambda (u) (- (pnhyper x (/ u)) alpha))
                            double-eps 1)))
                 (t 1))))
           (ncp-l (x alpha)
             (and (= x lo) 0)
             (let ((p (pnhyper x 1 :upper-tail t)))
               (cond 
                 ((> p alpha)
                  (funcall uniroot (lambda (u) (- (pnhyper x u :upper-tail t) alpha))
                         0 1))
                 ((< p alpha)
                  (/ (funcall uniroot (lambda (u) (- (pnhyper x (/ u) :upper-tail t) alpha))
                            double-eps 1)))
                 (t 1)))))
        (let ((p-value 
               (ecase alternative
                 (less (pnhyper x00 or))
                 (greater (pnhyper x00 or :upper-tail t)) 
                 (two-sided 
                  (let* ((relErr  (1+ 1.0d-7))
                         (d (dnhyper or))
                         (dstar (* (elt d (- x00 lo)) relErr)))
                    (loop for di in d when (< di dstar) summing di)))))
              (c-interval
               (if conf-int
                   (ecase alternative
                     (less (list 0 (ncp-u x00 (- 1 conf-level))))
                     (greater (list (ncp-l x00 (- 1 conf-level)) 'infinity))
                     (two-sided 
                      (let ((alpha (/ (- 1 conf-level) 2)))
                        (list (ncp-l x00 alpha) (ncp-u x00 alpha)))))
                   nil))
              (estimate (mle x00)))
          (values p-value c-interval estimate)))))
;;(fisher-test #2a((10 10) (10 20)))

lisp

; (DEFUN RIDDER ...) is being compiled.
;; The variable FNEW is undefined.
;; The compiler will assume this variable is a global.
Value = 17624

With output:

LatexWiki Image