|
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, Reduce, Maxima, Common Lisp and Sage (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 conditional maximum likelihood estimator for the odds ratio. The reference implementation is R, where the code and output would be:
As a caveat: I have little experience with these programs. Any changes or improvements to the programs would be welcomed.
To summarise, all six languages (Spad, Boot, Reduce, Maxima, Common Lisp and Sage) 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. Nested functions in Spad required the use of #1 and #2 argument references, although this has changed recently in Fricas; I appreciated Spad's lexical scoping and facility to fall back to a symbolic analysis. An Aldor implementation would be similar to 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 clumsy. 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. For Reduce, the lack of nested procedures in the algebraic mode made progress slow; importantly, the switch to symbolic mode made the implementation quite straightforward, with reasonably good debugging. The Maxima and Sage versions were particularly short, given that they already provide a root solver.
This begs the question: when would one use any of these languages for mixed numerical/symbolic analysis? In my opinion, 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. For an R user, Spad's type system seemed both fussy and extremely elegant; I found that debugging could be slow. Importantly, Lisp and Boot are able to evaluate Spad functions, potentially allowing the use of Lisp functions such as those translated from Fortran using f2cl (see also [SandboxMLE] for another example of numerical code in Spad). Axiom, Reduce, Maxima and Sage provide polished environments worthy of further consideration.
First, the Spad implementation:
fricas (1) -> <spad>
fricas )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
ridder(mnhyper(#1) - x::R,0,1)
else if mu<x::R then
1/ridder(mnhyper(1/#1) - x::R,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
ridder(pnhyper(x,#1,false) - alpha, 0.0::R, 1.0::R)
else if p>alpha then
1/ridder(pnhyper(x,1/#1,false) - alpha, 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
ridder(pnhyper(x,#1,true) - alpha, 0,1)
else if p<alpha then
1/ridder(pnhyper(x,1/#1,true) - alpha, 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>
fricas Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2045700826673834242-25px001.spad
using old system compiler.
TESTP abbreviates package TestPackage
------------------------------------------------------------------------
initializing NRLIB TESTP for TestPackage
compiling into NRLIB TESTP
importing TrigonometricFunctionCategory
compiling exported ridder : (Float -> Float,Float,Float) -> Float
Time: 0 SEC.
compiling exported msign : (Float,Float) -> Float
Time: 0 SEC.
compiling exported choose : (Integer,Integer) -> Fraction Integer
Time: 0 SEC.
compiling exported dhyper : (Integer,Integer,Integer,Integer) -> Fraction Integer
Time: 0 SEC.
compiling exported phyper : (Integer,Integer,Integer,Integer,Boolean) -> Fraction Integer
Time: 0 SEC.
compiling exported fisherTest : (Integer,Integer,Integer,Integer,String,Float,Boolean,Float) -> Record(PValue: Float,CI: List Float,Estimate: Float)
****** comp fails at level 8 with expression: ******
error in function fisherTest
(SEQ (|:=| (|:| |m| (|Integer|)) (+ |a| |c|))
(|:=| (|:| |n| (|Integer|)) (+ |b| |d|))
(|:=| (|:| |k| (|Integer|)) (+ |a| |b|))
(|:=| (|:| |x00| (|Integer|)) |a|)
(|:=| (|:| |lo| (|Integer|)) (|max| 0 (- |k| |n|)))
(|:=| (|:| |hi| (|Integer|)) (|min| |k| |m|))
(|:=| (|:| |support| (|List| (|Integer|)))
(COLLECT (IN |i| (SEGMENT |lo| |hi|)) |i|))
(|:=| (|:| |logdc| (|List| (|Float|)))
(COLLECT (IN |i| |support|)
(|log| (|::| (|dhyper| |i| |m| |n| |k|) (|Float|)))))
(|:=| (|:| |doubleEps| (|Float|))
(|::| ((|Sel| (|Float|) |float|) 1 -50 10) (|Float|)))
(|:=| (|:| |plusInfinity| (|Float|))
(|::| ((|Sel| (|Float|) |float|) 1 6 10) (|Float|)))
(|:=| (|:| |dnhyper| (|Mapping| (|List| (|Float|)) (|Float|)))
(SEQ (|:=| (|:| |ncp| (|Float|)) |#1|)
(|:=| (|:| |d| (|List| (|Float|)))
(COLLECT (IN |i| (SEGMENT 1 (|#| |logdc|)))
(+ (|logdc| |i|)
(* (|log| |ncp|) (|::| (|support| |i|) (|Float|))))))
(|:=| (|:| |maxd| (|Float|)) (|reduce| |max| |d|))
(|:=| (|:| |d2| (|List| (|Float|)))
(COLLECT (IN |di| |d|) (|exp| (- |di| |maxd|))))
(|:=| (|:| |sumd2| (|Float|)) (|reduce| "+" |d2|))
(|exit| 1 (COLLECT (IN |d2i| |d2|) (/ |d2i| |sumd2|)))))
(|:=| (|:| |mnhyper| (|Mapping| (|Float|) (|Float|)))
(SEQ (|:=| (|:| |ncp| (|Float|)) |#1|)
(|:=| (|:| #1=#:G12 (|Boolean|))
(= |ncp| (|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|))))
(|exit| 1
(IF #1#
(|::| |lo| (|Float|))
(SEQ (|:=| (|:| |d| (|List| (|Float|))) (|dnhyper| |ncp|))
(|exit| 1
(|reduce| "+"
(COLLECT (IN |i| (SEGMENT 1 (|#| |d|)))
(* (|::| (|support| |i|) (|Float|))
(|d| |i|))))))))))
(|:=|
(|:| |pnhyper| (|Mapping| (|Float|) (|Integer|) (|Float|) (|Boolean|)))
(SEQ (|:=| (|:| |q| (|Integer|)) |#1|)
(|:=| (|:| |ncp| (|Float|)) (|#| 2))
(|:=| (|:| |upperTail| (|Boolean|)) (|#| 3))
(|:=| (|:| #2=#:G14 (|Boolean|))
(= |ncp| ((|Sel| (|Float|) |float|) 1 0 10)))
(|exit| 1
(IF #2#
(IF |upperTail|
(|::| (|phyper| (- |q| 1) |m| |n| |k| |false|) (|Float|))
(|::| (|phyper| |q| |m| |n| |k| |true|) (|Float|)))
(SEQ
(|:=| (|:| #3=#:G13 (|Boolean|))
(= |ncp| ((|Sel| (|Float|) |float|) 0 0 10)))
(|exit| 1
(IF #3#
(IF |upperTail|
(IF (<= |q| |lo|)
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))
(|::| ((|Sel| (|Float|) |float|) 0 0 10)
(|Float|)))
(IF (>= |q| |lo|)
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))
(|::| ((|Sel| (|Float|) |float|) 0 0 10)
(|Float|))))
(SEQ
(|:=| (|:| |d| (|List| (|Float|))) (|dnhyper| |ncp|))
(|exit| 1
(IF |upperTail|
(|reduce| "+"
(COLLECT (IN |i| (SEGMENT 1 (|#| |d|)))
(|\|| (>= (|support| |i|) |q|))
(|d| |i|)))
(|reduce| "+"
(COLLECT (IN |i| (SEGMENT 1 (|#| |d|)))
(|\|| (<= (|support| |i|) |q|))
(|d| |i|)))))))))))))
(|:=| (|:| |mle| (|Mapping| (|Float|) (|Integer|)))
(SEQ (|:=| (|:| |x| (|Integer|)) |#1|)
(|exit| 1
(IF (= |x| |lo|)
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|))
(IF (= |x| |hi|)
|plusInfinity|
(SEQ
(|:=| (|:| |mu| (|Float|))
(|mnhyper|
(|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|))))
(|exit| 1
(IF (> |mu| (|::| |x| (|Float|)))
(|ridder| (- (|mnhyper| |#1|) (|::| |x| (|Float|))) 0
1)
(IF (< |mu| (|::| |x| (|Float|)))
(/ 1
(|ridder|
(- (|mnhyper| (/ 1 |#1|))
(|::| |x| (|Float|)))
|doubleEps|
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))))
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|)))))))))))
(|:=| (|:| |ncpU| (|Mapping| (|Float|) (|Integer|) (|Float|)))
(SEQ (|:=| (|:| |x| (|Integer|)) |#1|)
(|:=| (|:| |alpha| (|Float|)) (|#| 2))
(|exit| 1
(IF (= |x| |hi|)
|plusInfinity|
(SEQ
(|:=| (|:| |p| (|Float|))
(|pnhyper| |x|
(|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|))
|false|))
(|exit| 1
(IF (< |p| |alpha|)
(|ridder| (- (|pnhyper| |x| |#1| |false|) |alpha|)
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|))
(|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|)))
(IF (> |p| |alpha|)
(/ 1
(|ridder|
(- (|pnhyper| |x| (/ 1 |#1|) |false|) |alpha|)
|doubleEps|
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))))
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))))))))))
(|:=| (|:| |ncpL| (|Mapping| (|Float|) (|Integer|) (|Float|)))
(SEQ (|:=| (|:| |x| (|Integer|)) |#1|)
(|:=| (|:| |alpha| (|Float|)) (|#| 2))
(|exit| 1
(IF (= |x| |lo|)
(|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|))
(SEQ (|:=| (|:| |p| (|Float|)) (|pnhyper| |x| 1 |true|))
(|exit| 1
(IF (> |p| |alpha|)
(|ridder| (- (|pnhyper| |x| |#1| |true|) |alpha|) 0
1)
(IF (< |p| |alpha|)
(/ 1
(|ridder|
(- (|pnhyper| |x| (/ 1 |#1|) |true|) |alpha|)
|doubleEps|
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))))
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))))))))))
(|:=| (|:| |pValue| (|Float|))
(IF (= |alternative| "less")
(|pnhyper| |x00| OR |false|)
(IF (= |alternative| "greater")
(|pnhyper| |x00| OR |true|)
(IF (= |alternative| "two-sided")
(SEQ
(|:=| |relErr|
(+ 1 (|::| ((|Sel| (|Float|) |float|) 1 -7 10) (|Float|))))
(|:=| |dn| (|dnhyper| OR))
(|:=| |dstar| (* (|dn| (+ (- |x00| |lo|) 1)) |relErr|))
(|exit| 1
(|reduce| "+"
(COLLECT (IN |di| |dn|) (|\|| (< |di| |dstar|)) |di|))))
(- (|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|)))))))
(|:=| (|:| |cInterval| (|List| (|Float|)))
(IF |confInt|
(IF (= |alternative| "less")
(|construct| (|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|))
(|ncpU| |x00|
(-
(|::| ((|Sel| (|Float|) |float|) 1 0 10)
(|Float|))
|confLevel|)))
(IF (= |alternative| "greater")
(|construct|
(|ncpL| |x00|
(- (|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|))
|confLevel|))
|plusInfinity|)
(IF (= |alternative| "two-sided")
(SEQ (|:=| |alpha| (/ (- 1 |confLevel|) 2))
(|exit| 1
(|construct| (|ncpL| |x00| |alpha|)
(|ncpU| |x00| |alpha|))))
|noBranch|)))
(|construct| (- (|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|)))
(-
(|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|))))))
(|:=| |estimate| (|mle| |x00|))
(|exit| 1 (|construct| |pValue| |cInterval| |estimate|)))
****** level 8 ******
$x:= 2
$m:= (String)
$f:=
((((|ncp| #) (|q| # #) (|#3| #) (|#2| #) ...)))
>> Apparent user error:
no mode found for
#1Using this code in Axiom:
fricas -- test code is correct
alltests()
There are no library operations named alltests
Use HyperDoc Browse or issue
)what op alltests
to learn if there is any operation containing " alltests " in its
name.
Cannot find a no-argument definition or library operation named
alltests .
Second, 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
The file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1874789261577970310-25px003.boot
is needed but does not exist.Using the Boot code from Axiom:
fricas alltests()$Lisp
alltests is not a lisp function and so cannot be used with $Lisp.
Third, for Reduce (which was also my first Reduce program):
symbolic;
nil
on rounded;
nil | reduce |
symbolic procedure msign(x, y);
(abs x) * (if y>0.0 then 1.0 else if y<0.0 then -1.0 else 0.0);
msign | reduce |
%% Numerical root finding using Ridders method
%% (Exit criteria hacked: how can one return from a repeat .. until statement?)
symbolic procedure ridders(func, x1, x2);
begin scalar eps, maxit, fl, fh, xl, xh, ans, xnew, iterNum, fnew;
eps:= 1.0e-12;
maxit:= 100;
fl := funcall(func, x1);
fh := funcall(func, x2);
xl := x1;
xh := x2;
ans := -1.0e30;
xnew := 0.0e0;
iterNum := 0;
if (fl*fh) > 0.0 then rederr "Initial points are not either side of zero.";
if fl=0.0 then x1
else if fh=0.0 then x2
%if (fl*fh) < 0.0 then
else repeat begin scalar xm, fm, ss;
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:= funcall(func, ans);
%write(fnew);
%if fnew=0.0 then return ans;
if msign(fm,fnew) neq fm then begin;
xl:= xm;
fl:= fm;
xh:= ans;
fh:= fnew;
end
else if msign(fl, fnew) neq fl then begin;
xh:= ans;
fh:= fnew;
end
else if msign(fh, fnew) neq fh then begin;
xl:= ans;
fl:= fnew;
end;
iterNum:=iterNum+1;
if iterNum >=maxit then rederr "Maximum iterations exceeded";
%if verbose then write xl, xh, fl, fh;
end until abs(fnew)<eps or abs(xh-xl) <= eps;
return ans
end;
ridders | reduce |
symbolic procedure choose(n, x);
begin scalar total, denom;
total := 1.0;
for denom:=1:x do <<
total:=total/denom*(n-denom+1) >>;
return total
end;
choose | reduce |
symbolic procedure dhyper(x, m, n, k);
choose(m, x) * choose(n, k - x) / choose(m + n, k);
dhyper
procedure phyper(x, m, n, k, lowerTail);
if lowerTail then
for i:=1:x sum dhyper(i, m, n, k)
else
for i:=(x+1):k sum dhyper(i, m, n, k);
phyper | reduce |
symbolic procedure testTolerance(x, y, atol);
if abs(x-y) <= atol then t else nil;
testtolerance | reduce |
symbolic procedure test1(); testTolerance(2*ridders(function(cos),0,2),
cdr reval(algebraic pi),1e-8);
test1
symbolic procedure test2(); testTolerance(choose(100, 5), 75287520, 0);
test2
symbolic procedure test3(); testTolerance(dhyper(5, 10, 7, 8), 0.3628137, 1.0e-7);
test3
symbolic procedure test4(); testTolerance(log(dhyper(5, 10, 7, 8)),-1.013866, 1.0e-7);
test4
symbolic procedure test5(); testTolerance(phyper(5, 10, 7, 8, t),0.7821884, 1.0e-7);
test5
symbolic procedure test6(); testTolerance(phyper(5, 10, 7, 8, nil),0.2178116, 1.0e-7);
test6
symbolic procedure testSet1(); {test1(), test2(), test3(), test4(), test5(), test6()};
testset1 | reduce |
testSet1();
(t t t t t t) | reduce |
%%
symbolic procedure fisherTest(a,b,c,d,alternative,oddsratio,confLevel);
begin scalar m, n, k, x00, lo, hi, doubleEps,
dnhyper,mnhyper,pnhyper,mle,ncpu,ncpl,
support,d2,logdc,alpha,
pvalue,cinterval;
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 := for i:=lo:hi collect i;
logdc := for each i in support collect log(dhyper(i, m, n, k));
doubleEps := 1.0e-10;
dnhyper := function(lambda ncp;
begin scalar maxd, sumd2, d, d2;
d := for i:=1:length(logdc) collect nth(logdc,i)+log(ncp)*nth(support,i);
maxd := apply(function(max),d);
d2 :=for each di in d collect exp(di-maxd);
sumd2 := for each i in d2 sum i;
return for each d2i in d2 collect d2i/sumd2
end);
mnhyper := function(lambda ncp;
begin scalar d, value;
value := if ncp=0.0 then lo
%else if ncp equal 'plusInfinity then hi
else begin;
d := funcall(dnhyper,ncp);
return for i:=1:length(d) sum nth(support,i)*nth(d,i)
end;
return value
end);
pnhyper := function(lambda(q,ncp,upperTail);
if ncp=1.0 then
(if upperTail then phyper(q-1, m, n, k, nil)
else phyper(q, m, n, k, t))
else if ncp=0.0 then
(if upperTail then
(if q<=lo then 1.0 else 0.0)
else (if q>=lo then 1.0 else 0.0))
else if ncp equal 'plusInfinity then
(if upperTail then
(if q<=hi then 1.0 else 0.0)
else (if q>= hi then 1.0 else 0.0))
else
begin scalar d, value;
d := funcall(dnhyper,ncp);
value := if upperTail then
for i:=1:length(d) sum (if nth(support,i)>=q then nth(d,i) else 0)
else for i:=1:length(d) sum (if nth(support,i)<=q then nth(d,i) else 0);
return value
end);
ncpL := function(lambda alpha;
if x00=lo then 0
else begin scalar p,value,f1,f2;
f1 := function(lambda y; funcall(pnhyper,x00,y,t) - alpha);
f2 := function(lambda y; funcall(pnhyper,x00,1/y,t) - alpha);
p := funcall(pnhyper,x00, 1.0, t);
value := if p>alpha then
ridders(f1, doubleEps, 1) % zero bound caused problems
else if p<alpha then 1/ridders(f2, doubleEps, 1)
else 1;
return value;
end);
ncpU := function(lambda alpha;
if x00=hi then 'plusInfinity
else begin scalar p,value,f1,f2;
f1 := function(lambda y; funcall(pnhyper,x00,y,nil) - alpha);
f2 := function(lambda y; funcall(pnhyper,x00,1/y,nil) - alpha);
p := funcall(pnhyper,x00, 1.0, nil);
value := if p<alpha then ridders(f1, 0, 1)
else if p>alpha then 1/ridders(f2, doubleEps, 1)
else 1;
return value;
end);
pvalue :=
if alternative equal 'less then funcall(pnhyper,x00,oddsratio,nil)
else if alternative equal 'greater then funcall(pnhyper,x00,oddsratio,t)
else if alternative equal 'twosided then
begin scalar relerr,dstar,dn;
relErr:= 1+1.0e-7;
dn:= funcall(dnhyper,oddsratio);
dstar:= nth(dn,x00-lo+1)*relErr;
return for each di in dn sum (if di<dstar then di else 0)
end
else -1;
cInterval :=
if alternative equal 'less then
{0.0, funcall(ncpU,1-confLevel)}
else if alternative equal 'greater then
{funcall(ncpL,1-confLevel), 'plusInfinity}
else if alternative equal 'twosided then
{funcall(ncpL,(1-confLevel)/2), funcall(ncpU,(1-confLevel)/2)}
else {-1,-1}; % no match
mle :=
if x00=lo then 0.0
else if x00=hi then 'plusInfinity
else
begin scalar mu, value, f1, f2;
f1 := function(lambda y; funcall(mnhyper,y) - x00);
f2 := function(lambda y; funcall(mnhyper,1/y) - x00);
mu := funcall(mnhyper,1.0);
value := if mu>x00 then ridders(f1,0,1)
else if mu<x00 then
1/ridders(f2,doubleEps,1.0)
else 1.0;
return value;
end;
return {pvalue,cinterval,mle};
end;
fishertest | reduce |
fishertest(10,10,10,20,'twosided,1,0.95);
(0.257549 (0.538399 7.43634) 1.97163)
fishertest(10,10,10,20,'less,1,0.95);
(0.929481 (0.0 6.14381) 1.97163)
fishertest(10,10,10,20,'greater,1,0.95);
(0.188301 (0.645994 plusinfinity) 1.97163) | reduce |
symbolic procedure testSet2;
begin scalar fit1;
fit1 := fishertest(10,10,10,20,'twosided,1,0.95);
return {testTolerance(first(fit1), 0.2575, 1.0e-3),
testTolerance(caadr(fit1),0.5383996, 1.0e-6),
testTolerance(cadadr(fit1),7.4363242, 1.0e-4),
testTolerance(third(fit1),1.971640, 1.0e-4)}
end;
testset2 | reduce |
testSet2();
(t t t t) | reduce |
Fourth, the Maxima implementation was fairly brief:
Fifth, 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 Your user access level is compiler and this command is therefore not
available. See the )set userlevel command for more information.With output:
Finally, I have also included an implementation in Sage and commented out early code that would
allow this to be used in Python.
![\begin{sageblock}
import scipy
from scipy.optimize import brentq</p>
<p>## the following is required for use in Python/Scipy outside of Sage
# from math import log, exp, cos, pi
# infinity="infinity"
# def binomial(n,x):
# total=1
# for i in range(min(x,n-x)):
# total=total*(n-i)/(i+1)
# return total</p>
<p>def dhyper(x, m, n, k):
return 1.0 <em> binomial(m, x) </em> binomial(n, k - x) / binomial(m + n, k)
def phyper(x, m, n, k, lowerTail):
if lowerTail:
return sum([dhyper(i, m, n, k) for i in range(1,x+1)]<a class=]( "create this page") ?)
else:
return sum([dhyper(i, m, n, k) for i in range(x+1,k+1)]?)
def fishertest(a,b,c,d,alternative,oddsratio,confLevel):
m = a+c
n = b+d
k = a+b
x00 = a
lo = max(0, k-n)
hi = min(k, m)
support = range(lo,hi+1)
logdc = [log(dhyper(i, m, n, k)) for i in support]?
doubleEps = 1.0e-10
def dnhyper(ncp):
d = [logdc[i]?+log(ncp)*support[i]? for i in range(len(logdc))]
d2 = [exp(di-max(d)) for di in d]?
return [d2i/sum(d2) for d2i in d2]?
def mnhyper(ncp):
if ncp==0:
return lo
if ncp==infinity:
return hi
else:
d = dnhyper(ncp)
return sum([support[i]?*d[i]? for i in range(len(d))])
def pnhyper(q,ncp,upperTail):
if ncp==1:
if upperTail:
return phyper(q-1, m, n, k, False)
else:
return phyper(q, m, n, k, True)
elif ncp==0:
if upperTail:
return 1 if q<=lo else 0
else:
return 1 if q>=lo else 0
elif ncp==infinity:
if upperTail:
return 1 if q<=hi else 0
else:
return 1 if q>= hi else 0
else:
d = dnhyper(ncp)
if upperTail:
return sum([d[i]? for i in range(len(d)) if support[i]?>=q])
else:
return sum([d[i]? for i in range(1,len(d)) if support[i]?<=q])
def ncpL(alpha):
if x00==lo:
return 0
else:
p = pnhyper(x00, 1, True)
if p>alpha:
return brentq(lambda y: pnhyper(x00,y,True) - alpha, 0, 1)
elif p<alpha:
return 1/brentq(lambda y: pnhyper(x00,1/y,True) - alpha, doubleEps, 1)
else:
return 1
def ncpU(alpha):
if x00==hi:
return infinity
else:
p = pnhyper(x00, 1, False)
if p<alpha:
return brentq(lambda y: pnhyper(x00,y,False) - alpha, 0, 1)
elif p>alpha:
return 1/brentq(lambda y: pnhyper(x00,1/y,False) - alpha, doubleEps, 1)
else:
return 1
def pvalue():
if alternative == "less":
return pnhyper(x00,oddsratio,False)
elif alternative == "greater":
return pnhyper(x00,oddsratio,True)
elif alternative == "two-sided":
relErr=1+1.0e-7
dn=dnhyper(oddsratio)
dstar=dn[x00-lo]?*relErr
return sum([di for di in dn if di<dstar]?)
else:
return -1
def cInterval():
if alternative == "less":
return [0, ncpU(1-confLevel)]?
elif alternative == "greater":
return [ncpL(1-confLevel), infinity]?
elif alternative == "two-sided":
return [ncpL((1-confLevel)/2.0), ncpU((1-confLevel)/2.0)]?
else:
return [-1,-1]?
def mle():
if x00==lo:
return 0
elif x00==hi:
return infinity
else:
mu = mnhyper(1)
if mu>x00:
return brentq(lambda y: mnhyper(y) - x00,0,1)
elif mu<x00:
return 1/brentq(lambda y: mnhyper(1/y) - x00,doubleEps,1)
else:
return 1
return dict(pvalue=pvalue(),cInterval=cInterval(),mle=mle())
## do some checks
def testTolerance(x,y,tol):
return abs(x-y)<=tol
fit=fishertest(10,10,10,20,"two-sided",1,0.95)
print [testTolerance(brentq(cos,0,2)*2,pi,0),
testTolerance(binomial(100, 5), 75287520, 0),
testTolerance(binomial(100, 95), 75287520, 0),
testTolerance(dhyper(5, 10, 7, 8), 0.3628137, 1.0e-7),
testTolerance(log(dhyper(5, 10, 7, 8)),-1.013866, 1.0e-7),
testTolerance(phyper(5, 10, 7, 8, True),0.7821884, 1.0e-7),
testTolerance(phyper(5, 10, 7, 8, False),0.2178116, 1.0e-7),
testTolerance(fit['pvalue']?, 0.2575, 1.0e-3),
testTolerance(fit['cInterval']?[0]?, 0.5383996, 1.0e-6),
testTolerance(fit['cInterval']?[1]?, 7.4363242, 1.0e-4),
testTolerance(fit['mle']?, 1.971640, 1.0e-4)]
## print some tests
print fishertest(10,10,10,20,"two-sided",1,0.95)
print fishertest(10,10,10,20,"less",1,0.95)
print fishertest(10,10,10,20,"greater",1,0.95)
" title="
\begin{sageblock}
import scipy
from scipy.optimize import brentq
## the following is required for use in Python/Scipy outside of Sage
# from math import log, exp, cos, pi
# infinity="infinity"
# def binomial(n,x):
# total=1
# for i in range(min(x,n-x)):
# total=total*(n-i)/(i+1)
# return total
def dhyper(x, m, n, k):
return 1.0 binomial(m, x) binomial(n, k - x) / binomial(m + n, k)
def phyper(x, m, n, k, lowerTail):
if lowerTail:
return sum([dhyper(i, m, n, k) for i in range(1,x+1)]?)
else:
return sum([dhyper(i, m, n, k) for i in range(x+1,k+1)]?)
def fishertest(a,b,c,d,alternative,oddsratio,confLevel):
m = a+c
n = b+d
k = a+b
x00 = a
lo = max(0, k-n)
hi = min(k, m)
support = range(lo,hi+1)
logdc = [log(dhyper(i, m, n, k)) for i in support]?
doubleEps = 1.0e-10
def dnhyper(ncp):
d = [logdc[i]?+log(ncp)*support[i]? for i in range(len(logdc))]
d2 = [exp(di-max(d)) for di in d]?
return [d2i/sum(d2) for d2i in d2]?
def mnhyper(ncp):
if ncp==0:
return lo
if ncp==infinity:
return hi
else:
d = dnhyper(ncp)
return sum([support[i]?*d[i]? for i in range(len(d))])
def pnhyper(q,ncp,upperTail):
if ncp==1:
if upperTail:
return phyper(q-1, m, n, k, False)
else:
return phyper(q, m, n, k, True)
elif ncp==0:
if upperTail:
return 1 if q<=lo else 0
else:
return 1 if q>=lo else 0
elif ncp==infinity:
if upperTail:
return 1 if q<=hi else 0
else:
return 1 if q>= hi else 0
else:
d = dnhyper(ncp)
if upperTail:
return sum([d[i]? for i in range(len(d)) if support[i]?>=q])
else:
return sum([d[i]? for i in range(1,len(d)) if support[i]?<=q])
def ncpL(alpha):
if x00==lo:
return 0
else:
p = pnhyper(x00, 1, True)
if p>alpha:
return brentq(lambda y: pnhyper(x00,y,True) - alpha, 0, 1)
elif p<alpha:
return 1/brentq(lambda y: pnhyper(x00,1/y,True) - alpha, doubleEps, 1)
else:
return 1
def ncpU(alpha):
if x00==hi:
return infinity
else:
p = pnhyper(x00, 1, False)
if p<alpha:
return brentq(lambda y: pnhyper(x00,y,False) - alpha, 0, 1)
elif p>alpha:
return 1/brentq(lambda y: pnhyper(x00,1/y,False) - alpha, doubleEps, 1)
else:
return 1
def pvalue():
if alternative == "less":
return pnhyper(x00,oddsratio,False)
elif alternative == "greater":
return pnhyper(x00,oddsratio,True)
elif alternative == "two-sided":
relErr=1+1.0e-7
dn=dnhyper(oddsratio)
dstar=dn[x00-lo]?*relErr
return sum([di for di in dn if di<dstar]?)
else:
return -1
def cInterval():
if alternative == "less":
return [0, ncpU(1-confLevel)]?
elif alternative == "greater":
return [ncpL(1-confLevel), infinity]?
elif alternative == "two-sided":
return [ncpL((1-confLevel)/2.0), ncpU((1-confLevel)/2.0)]?
else:
return [-1,-1]?
def mle():
if x00==lo:
return 0
elif x00==hi:
return infinity
else:
mu = mnhyper(1)
if mu>x00:
return brentq(lambda y: mnhyper(y) - x00,0,1)
elif mu<x00:
return 1/brentq(lambda y: mnhyper(1/y) - x00,doubleEps,1)
else:
return 1
return dict(pvalue=pvalue(),cInterval=cInterval(),mle=mle())
## do some checks
def testTolerance(x,y,tol):
return abs(x-y)<=tol
fit=fishertest(10,10,10,20,"two-sided",1,0.95)
print [testTolerance(brentq(cos,0,2)*2,pi,0),
testTolerance(binomial(100, 5), 75287520, 0),
testTolerance(binomial(100, 95), 75287520, 0),
testTolerance(dhyper(5, 10, 7, 8), 0.3628137, 1.0e-7),
testTolerance(log(dhyper(5, 10, 7, 8)),-1.013866, 1.0e-7),
testTolerance(phyper(5, 10, 7, 8, True),0.7821884, 1.0e-7),
testTolerance(phyper(5, 10, 7, 8, False),0.2178116, 1.0e-7),
testTolerance(fit['pvalue']?, 0.2575, 1.0e-3),
testTolerance(fit['cInterval']?[0]?, 0.5383996, 1.0e-6),
testTolerance(fit['cInterval']?[1]?, 7.4363242, 1.0e-4),
testTolerance(fit['mle']?, 1.971640, 1.0e-4)]
## print some tests
print fishertest(10,10,10,20,"two-sided",1,0.95)
print fishertest(10,10,10,20,"less",1,0.95)
print fishertest(10,10,10,20,"greater",1,0.95)
" class="equation" src="images/7581108121597389083-16.0px.png" align="bottom" Style="vertical-align:text-bottom" width="816" height="1056"/>
This would give the following output:
[0 <= 0, True, True, True, True, True, True, 4.92428109829e-05 <= 0.00100000000000000,
True, True, True]
{'mle': 1.9716269432603342, 'pvalue': 0.257549242811,
'cInterval': [0.53839938167817269, 7.4363408387619616]}
{'mle': 1.9716269432603342, 'pvalue': 0.929480494158038,
'cInterval': [0, 6.1438085831669706]}
{'mle': 1.9716269432603342, 'pvalue': 0.188301375769915,
'cInterval': [0.64599378194777934, +Infinity]}
Some or all expressions may not have rendered properly,
because Maxima returned the following error:
Maxima error: Killed
ulimit -t 240; /usr/local/bin/maxima -p /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/../../Products/ZWiki/plugins/mathaction/mathaction-maxima-5.30.0.lisp < /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/5516975315810540791-25px.mbat
<maxima>
Maxima 5.30.0 http://maxima.sourceforge.net
using Lisp SBCL 2.2.9.debian
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
<prompt>(\mathrm{\%i}1) </prompt>
read and interpret file: /var/aw/var/LatexWiki/2074375897756092749-25px.001.max
(%i2) load(distrib)
<latex>\mbox{\verb|f|}\mbox{\verb|i|}\mbox{\verb|l|}\mbox{\verb|e|}
\mbox{\verb|_|}\mbox{\verb|s|}\mbox{\verb|e|}\mbox{\verb|a|}
\mbox{\verb|r|}\mbox{\verb|c|}\mbox{\verb|h|}\mbox{\verb|1|}
\mbox{\verb|:|}\mbox{\verb| |}\mathrm{distrib}\mbox{\verb| |}
\mbox{\verb|n|}\mbox{\verb|o|}\mbox{\verb|t|}\mbox{\verb| |}
\mbox{\verb|f|}\mbox{\verb|o|}\mbox{\verb|u|}\mbox{\verb|n|}
\mbox{\verb|d|}\mbox{\verb| |}\mbox{\verb|i|}\mbox{\verb|n|}
\mbox{\verb| |}\mbox{\verb|f|}\mbox{\verb|i|}\mbox{\verb|l|}
\mbox{\verb|e|}\mbox{\verb|_|}\mbox{\verb|s|}\mbox{\verb|e|}
\mbox{\verb|a|}\mbox{\verb|r|}\mbox{\verb|c|}\mbox{\verb|h|}
\mbox{\verb|_|}\mbox{\verb|m|}\mbox{\verb|a|}\mbox{\verb|x|}
\mbox{\verb|i|}\mbox{\verb|m|}\mbox{\verb|a|}\mbox{\verb|,|}
\mbox{\verb|f|}\mbox{\verb|i|}\mbox{\verb|l|}\mbox{\verb|e|}
\mbox{\verb|_|}\mbox{\verb|s|}\mbox{\verb|e|}\mbox{\verb|a|}
\mbox{\verb|r|}\mbox{\verb|c|}\mbox{\verb|h|}\mbox{\verb|_|}
\mbox{\verb|l|}\mbox{\verb|i|}\mbox{\verb|s|}\mbox{\verb|p|}
\mbox{\verb|.|}</latex>
-- an error. To debug this try: debugmode(true);
<latex>\mbox{\tt\red(\mathrm{\%o2}) \black}
\mbox{{}/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2074375897756092749-25px.001.max{}}
</latex>
<prompt>(\mathrm{\%i}3) </prompt>
read and interpret file: /var/aw/var/LatexWiki/4805595149220346197-25px.002.max
(%i4) numer:true
<latex>\mbox{\tt\red(\mathrm{\%o4}) \black}\mathbf{true}
</latex><latex>\mbox{\tt\red(\mathrm{\%o4}) \black}
\mbox{{}/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/4805595149220346197-25px.002.max{}}
</latex>
<prompt>(\mathrm{\%i}5) </prompt>
read and interpret file: /var/aw/var/LatexWiki/2094584513860488517-25px.003.max
(%i6) display2d:false
(%i7) fishertest(a,b,c,d,alternative,oddsratio,confLevel):=block(
[m,n,k,x00,lo,hi,doubleEps,dnhyper,mnhyper,pnhyper,mle,ncpu,
ncpl,support,d2,logdc,alpha,pvalue,cinterval],m:c+a,n:d+b,
k:b+a,x00:a,lo:max(0,k-n),hi:min(k,m),
support:makelist(i,i,lo,hi),
logdc:create_list(log(pdf_hypergeometric(i,m,n,k)),i,
support),doubleEps:1.0e-10,
dnhyper:lambda([ncp],
block([maxd,sumd2,d,d2],
d:makelist(log(ncp)*support[i]+logdc[i],
i,1,length(logdc)),
maxd:apply(max,d),
d2:create_list(exp(di-maxd),di,d),
sumd2:lsum(i,i,d2),
create_list(d2i/sumd2,d2i,d2))),
mnhyper:lambda([ncp],
block([d,value],
if equal(ncp,0) then lo
else (if equal(ncp,inf) then hi
else block([],
d:dnhyper(ncp),
sum(
support[i]
*d[i],i,1,
length(
d)))))),
pnhyper:lambda([q,ncp,upperTail],
if equal(ncp,1)
then (if upperTail
then 1
-cdf_hypergeometric(q-1,m,n,k)
else cdf_hypergeometric(
q,m,n,k))
else (if equal(ncp,0)
then (if upperTail
then (if q <= lo
then 1
else 0)
else (if q >= lo
then 1
else 0))
else (if equal(ncp,inf)
then (if upperTail
then (if
q <= hi
then 1
else 0)
else (if
q >= hi
then 1
else 0))
else block(
[d],d:dnhyper(ncp),
if upperTail
then sum(
if support[i] >= q
then d[i]
else 0,i,1,
length(d))
else sum(
if support[i] <= q
then d[i]
else 0,i,1,
length(d)))))),
ncpL:lambda([alpha],
if equal(x00,lo) then 0
else block([p],p:pnhyper(x00,1,true),
if p > alpha
then find_root(
lambda([y],
pnhyper(x00,y,true)
-alpha),0,1)
else (if p < alpha
then 1
/find_root(
lambda([y],
pnhyper(
x00,1/y,
true)
-alpha),
doubleEps,1)
else 1))),
ncpU:lambda([alpha],
if x00 = hi then inf
else block([p],p:pnhyper(x00,1,false),
if p < alpha
then find_root(
lambda([y],
pnhyper(x00,y,false)
-alpha),0,1)
else (if p > alpha
then 1
/find_root(
lambda([y],
pnhyper(
x00,1/y,
false)
-alpha),
doubleEps,1)
else 1))),
pvalue:if alternative = 'less
then pnhyper(x00,oddsratio,false)
else (if alternative = 'greater
then pnhyper(x00,oddsratio,true)
else (if alternative = 'twosided
then block([relerr,dstar,dn],
relErr:1.0e-7+1,
dn
:dnhyper(
oddsratio),
dstar
:dn[1-lo+x00]
*relErr,
lsum(
if di < dstar
then di
else 0,di,dn))
else -1)),
cInterval:if alternative = 'less then [0,ncpU(1-confLevel)]
else (if alternative = 'greater
then [ncpL(1-confLevel),inf]
else (if alternative = 'twosided
then [
ncpL((1-confLevel)/2),
ncpU((1-confLevel)/2)]
else [-1,-1])),
mle:if equal(x00,lo) then 0
else (if equal(x00,hi) then inf
else block([mu],mu:mnhyper(1),
if mu > x00
then find_root(
lambda([y],mnhyper(y)-x00),
0,1)
else (if mu < x00
then 1
/find_root(
lambda([y],
mnhyper(
1/y)
-x00),
doubleEps,1)
else 1))),
[pvalue,cInterval,mle])
(%i8) display2d:true
<latex>\mbox{\tt\red(\mathrm{\%o8}) \black}
\mbox{{}/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/2094584513860488517-25px.003.max{}}
</latex>
<prompt>(\mathrm{\%i}9) </prompt>
read and interpret file: /var/aw/var/LatexWiki/8899926541800781278-25px.004.max
(%i10) fishertest(10,10,10,20,'twosided,1,0.95)
|