|
|
|
last edited 13 years ago by Mark Clements |
Edit detail for SandBoxFisher revision 6 of 11
| 1 2 3 4 5 6 7 8 9 10 11 | ||
|
Editor: Mark Clements
Time: 2009/04/07 01:14:05 GMT-7 |
||
| Note: | ||
added:
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 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}
As a caveat: I have little experience with these programs. Any changes or improvements to the programs would be welcomed.
To summarise, all five languages (Spad, Boot, Reduce, Maxima 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 clumsy, with the requirement to use the #1 and #2 argument references (although this has changed recently in Fricas); I did, however, appreciate Spad's lexical scoping and facility to fall back to a symbolic analysis. An Aldor implementation may 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 particularly 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 version was particularly short, given that Maxima already provides a root solver and the hypergeometric distribution.
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 Reduce and 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 extremely elegant); I also found that debugging could be slow. The lack of ability to call Spad functions from Lisp or Boot is an important 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). Reduce and Maxima both provided fairly polished environments worthy of further consideration.
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 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> fricasCompiling FriCAS source code from file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/303933172962875058-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.01 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|))) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|mnhyper| |#1|) (|::| |x| (|Float|)))) (|exit| 1 (|ridder| |f| 0 1))) (IF (< |mu| (|::| |x| (|Float|))) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|mnhyper| (/ 1 |#1|)) (|::| |x| (|Float|)))) (|exit| 1 (/ 1 (|ridder| |f| |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|) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|pnhyper| |x| |#1| |false|) |alpha|)) (|exit| 1 (|ridder| |f| (|::| ((|Sel| (|Float|) |float|) 0 0 10) (|Float|)) (|::| ((|Sel| (|Float|) |float|) 1 0 10) (|Float|))))) (IF (> |p| |alpha|) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|pnhyper| |x| (/ 1 |#1|) |false|) |alpha|)) (|exit| 1 (/ 1 (|ridder| |f| |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|) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|pnhyper| |x| |#1| |true|) |alpha|)) (|exit| 1 (|ridder| |f| 0 1))) (IF (< |p| |alpha|) (SEQ (|:=| (|:| |f| (|Mapping| (|Float|) (|Float|))) (- (|pnhyper| |x| (/ 1 |#1|) |true|) |alpha|)) (|exit| 1 (/ 1 (|ridder| |f| |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 .The Boot translation was more fiddly - but, then again, I had never used Boot before.
bootdoubleFloat(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()] bootThe file /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1874789261577970310-25px003.boot is needed but does not exist.Using this code in Axiom:
fricasalltests()$Lisp
alltests is not a lisp function and so cannot be used with $Lisp.For Reduce (which was also my first Reduce program):
symbolic; | 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); | 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; | 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; | reduce |
symbolic procedure dhyper(x, m, n, k); choose(m, x) * choose(n, k - x) / choose(m + n, k); | reduce |
symbolic procedure testTolerance(x, y, atol); if abs(x-y) <= atol then t else nil; | reduce |
symbolic procedure test1(); testTolerance(2*ridders(function(cos),0,2),
cdr reval(algebraic pi),1e-8); | reduce |
testSet1(); | 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; | reduce |
fishertest(10,10,10,20,'twosided,1,0.95); | 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; | reduce |
testSet2(); | reduce |
The Maxima implementation was fairly brief:
\begin{maxima}
load(distrib);
\end{maxima}
\begin{maxima}
numer:true;
\end{maxima}
\begin{maxima}
display2d:false$
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 : a+c,
n : b+d,
k : a+b,
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(logdc[i]?+log(ncp)*support[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 The implementation in Common Lisp was a more direct translation of the R code:;; from cl-statistics.lisp
(defun safe-exp (x)
"Eliminates floating point underflow for the exponential function.
Instead,
(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))) Your user access level is compiler and this command is therefore not
available. See the )set userlevel command for more information.
With output:
\begin{verbatim} CL-USER> (fisher-test #2a((10 10) (10 20))) 0.2575492428109829d0 (0.5383993816781727d0 7.4363408387439875d0) 1.9716269432603386d0 \end{verbatim}
Some or all expressions may not have rendered properly, because Maxima returned the following error:
Maxima error: Killedulimit -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)
Some or all expressions may not have rendered properly, because Latex returned the following error:
This is pdfTeX, Version 3.141592653-2.6-1.40.24 (TeX Live 2022/Debian) (preloaded format=latex) restricted \write18 enabled. entering extended mode (./8783698267322362323-16.0px.tex LaTeX2e <2022-11-01> patch level 1 L3 programming layer <2023-01-16> (/usr/share/texlive/texmf-dist/tex/latex/base/article.cls Document Class: article 2022/07/02 v1.4n Standard LaTeX document class (/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo)) (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucs.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/data/uni-global.def)) (/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty (/usr/share/texlive/texmf-dist/tex/latex/ucs/utf8x.def)) (/usr/share/texlive/texmf-dist/tex/latex/bbm-macros/bbm.sty) (/usr/share/texlive/texmf-dist/tex/latex/jknapltx/mathrsfs.sty) (/usr/share/texlive/texmf-dist/tex/latex/base/fontenc.sty) (/usr/share/texlive/texmf-dist/tex/latex/pstricks/pstricks.sty (/usr/share/texlive/texmf-dist/tex/generic/iftex/iftex.sty) (/usr/share/texlive/texmf-dist/tex/latex/xcolor/xcolor.sty (/usr/share/texlive/texmf-dist/tex/latex/graphics-cfg/color.cfg) (/usr/share/texlive/texmf-dist/tex/latex/graphics-def/dvips.def) (/usr/share/texlive/texmf-dist/tex/latex/graphics/mathcolor.ltx) (/usr/share/texlive/texmf-dist/tex/latex/graphics/dvipsnam.def)) (/usr/share/texlive/texmf-dist/tex/generic/xkeyval/pst-xkey.tex (/usr/share/texlive/texmf-dist/tex/latex/xkeyval/xkeyval.sty (/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkeyval.tex (/usr/share/texlive/texmf-dist/tex/generic/xkeyval/xkvutils.tex (/usr/share/texlive/texmf-dist/tex/generic/xkeyval/keyval.tex))))) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks.tex (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pst-fp.tex `pst-fp' v0.06, 2020/11/20 (hv)) (/usr/share/texlive/texmf-dist/tex/latex/pgf/utilities/pgffor.sty (/usr/share/texlive/texmf-dist/tex/latex/pgf/utilities/pgfrcs.sty (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgfutil-common.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgfutil-latex.def) (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgfrcs.code.tex (/usr/share/texlive/texmf-dist/tex/generic/pgf/pgf.revision.tex))) (/usr/share/texlive/texmf-dist/tex/latex/pgf/utilities/pgfkeys.sty (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgfkeys.code.tex (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgfkeyslibraryfiltered .code.tex))) (/usr/share/texlive/texmf-dist/tex/latex/pgf/math/pgfmath.sty (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathutil.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathparser.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.basic.code .tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.trigonomet ric.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.random.cod e.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.comparison .code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.base.code. tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.round.code .tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.misc.code. tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.integerari thmetics.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathcalc.code.tex) (/usr/share/texlive/texmf-dist/tex/generic/pgf/math/pgfmathfloat.code.tex))) (/usr/share/texlive/texmf-dist/tex/generic/pgf/utilities/pgffor.code.tex)) `PSTricks' v3.18 <2022/11/28> (tvz,hv) --- We are running latex or xelatex --- (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks.con) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks-color.tex) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks-arrows.tex) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks-dots.tex) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pstricks.con)) (/usr/share/texlive/texmf-dist/tex/generic/pstricks/pst-fp.tex `pst-fp' v0.06, 2020/11/20 (hv))) (/usr/share/texlive/texmf-dist/tex/latex/graphics/epsfig.sty (/usr/share/texlive/texmf-dist/tex/latex/graphics/graphicx.sty (/usr/share/texlive/texmf-dist/tex/latex/graphics/graphics.sty (/usr/share/texlive/texmf-dist/tex/latex/graphics/trig.sty) (/usr/share/texlive/texmf-dist/tex/latex/graphics-cfg/graphics.cfg)))) (/usr/share/texlive/texmf-dist/tex/latex/pst-grad/pst-grad.sty (/usr/share/texlive/texmf-dist/tex/generic/pst-grad/pst-grad.tex `pst-grad' v1.06, 2006/11/27 (tvz,dg,hv))) (/usr/share/texlive/texmf-dist/tex/latex/pst-plot/pst-plot.sty (/usr/share/texlive/texmf-dist/tex/latex/xkeyval/pst-xkey.sty) (/usr/share/texlive/texmf-dist/tex/latex/multido/multido.sty (/usr/share/texlive/texmf-dist/tex/generic/multido/multido.tex v1.42, 2010/05/14 <tvz>)) (/usr/share/texlive/texmf-dist/tex/generic/pst-plot/pst-plot.tex (/usr/share/texlive/texmf-dist/tex/generic/pst-tools/pst-tools.tex `PST-tools' v0.12, 2021/09/23 (hv)) (/usr/share/texlive/texmf-dist/tex/generic/pstricks-add/pstricks-add.tex (/usr/share/texlive/texmf-dist/tex/generic/pst-node/pst-node.tex v1.43, 2022/01/31) (/usr/share/texlive/texmf-dist/tex/generic/pst-arrow/pst-arrow.tex `pst-arrow' v0.05, 2021/11/16 (dr,hv)) (/usr/share/texlive/texmf-dist/tex/generic/pst-3d/pst-3d.tex `PST-3d' v1.11, 2010/02/14 (tvz)) (/usr/share/texlive/texmf-dist/tex/generic/pst-math/pst-math.tex `pst-math' v0.66 , (CJ,hv)) `pstricks-add' v3.93, 2022/11/21 (dr,hv)) v1.94, 2022/11/21 (tvz,hv))) (/usr/share/texlive/texmf-dist/tex/latex/geometry/geometry.sty (/usr/share/texlive/texmf-dist/tex/generic/iftex/ifvtex.sty)Package geometry Warning: `lmargin' and `rmargin' result in NEGATIVE (-108.405p t). `width' should be shortened in length.
) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsmath.sty For additional information on amsmath, use the `?' option. (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amstext.sty (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsgen.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsbsy.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsmath/amsopn.sty)) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amsfonts.sty) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/amssymb.sty) (/usr/share/texlive/texmf-dist/tex/latex/amscls/amsthm.sty) (/usr/share/texlive/texmf-dist/tex/latex/setspace/setspace.sty) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xy.sty (/usr/share/texlive/texmf-dist/tex/generic/xypic/xy.tex Bootstrap'ing: catcodes, docmode, (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyrecat.tex ) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyidioms.tex)
Xy-pic version 3.8.9 <2013/10/06> Copyright (c) 1991-2013 by Kristoffer H. Rose <krisrose@tug.org> and others Xy-pic is free software: see the User's Guide for details.
Loading kernel: messages; fonts; allocations: state, direction, utility macros; pictures: \xy, positions, objects, decorations; kernel objects: directionals, circles, text; options; algorithms: directions, edges, connections; Xy-pic loaded) (/usr/share/texlive/texmf-dist/tex/generic/iftex/ifpdf.sty)) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyall.tex Xy-pic option: All features v.3.8 (/usr/share/texlive/texmf-dist/tex/generic/xypic/xycurve.tex Xy-pic option: Curve and Spline extension v.3.12 curve, circles, loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyframe.tex Xy-pic option: Frame and Bracket extension v.3.14 loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xycmtip.tex Xy-pic option: Computer Modern tip extension v.3.7 (/usr/share/texlive/texmf-dist/tex/generic/xypic/xytips.tex Xy-pic option: More Tips extension v.3.11 loaded) loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyline.tex Xy-pic option: Line styles extension v.3.10 loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyrotate.tex Xy-pic option: Rotate and Scale extension v.3.8 loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xycolor.tex Xy-pic option: Colour extension v.3.11 loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xymatrix.tex Xy-pic option: Matrix feature v.3.14 loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xyarrow.tex Xy-pic option: Arrow and Path feature v.3.9 path, \ar, loaded) (/usr/share/texlive/texmf-dist/tex/generic/xypic/xygraph.tex Xy-pic option: Graph feature v.3.11 loaded) loaded) (/usr/share/texlive/texmf-dist/tex/latex/tools/verbatim.sty) (/usr/share/texlive/texmf-dist/tex/latex/graphviz/graphviz.sty (/usr/share/texlive/texmf-dist/tex/latex/psfrag/psfrag.sty)) (/usr/share/texmf/tex/latex/sagetex.sty Writing sage input file 8783698267322362323-16.0px.sage ) (/usr/share/texlive/texmf-dist/tex/latex/gnuplottex/gnuplottex.sty (/usr/share/texlive/texmf-dist/tex/latex/moreverb/moreverb.sty) (/usr/share/texlive/texmf-dist/tex/latex/base/ifthen.sty) (/usr/share/texlive/texmf-dist/tex/generic/catchfile/catchfile.sty (/usr/share/texlive/texmf-dist/tex/generic/infwarerr/infwarerr.sty) (/usr/share/texlive/texmf-dist/tex/generic/ltxcmds/ltxcmds.sty) (/usr/share/texlive/texmf-dist/tex/generic/etexcmds/etexcmds.sty))
Package gnuplottex Warning: Shell escape not enabled. (gnuplottex) You'll need to convert the graphs yourself.
) (/usr/share/texlive/texmf-dist/tex/latex/l3backend/l3backend-dvips.def) No file 8783698267322362323-16.0px.aux. (/usr/share/texlive/texmf-dist/tex/latex/ucs/ucsencs.def) geometry driver: auto-detecting geometry detected driver: dvips (/usr/share/texlive/texmf-dist/tex/latex/base/t1cmtt.fd)
LaTeX Warning: Characters dropped after `\end{verbatim}' on input line 135.
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.136 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.138 \end{maxima} \newpage [1]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.139 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.141 \end{maxima} \newpage [2]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.142 \begin{maxima}
(/usr/share/texlive/texmf-dist/tex/latex/jknapltx/ursfs.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsa.fd) (/usr/share/texlive/texmf-dist/tex/latex/amsfonts/umsb.fd)
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.239 \end{maxima} \newpage Missing $ inserted. <inserted text> $ l.239 \end{maxima} \newpage
Overfull \hbox (321.19646pt too wide) in paragraph at lines 143--239 []\T1/cmr/m/n/12 display2d:false$\OML/cmm/m/it/12 fishertest\OT1/cmr/m/n/12 (\O ML/cmm/m/it/12 a; b; c; d; alternative; oddsratio; confLevel\OT1/cmr/m/n/12 ) : = \OML/cmm/m/it/12 block\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 m; n; k; x\OT1/cmr/m /n/12 00\OML/cmm/m/it/12 ; lo; hi; doubleEps; dnhyper; mnhyper; pnhyper; mle; n cpu; ncpl; support; d\OT1/cmr/m/n/12 2\OML/cmm/m/it/12 ; logdc; alpha; pvalue; cinterval\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; m \OT1/cmr/m/n/12 :
Overfull \hbox (137.83548pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 a \OT1/cmr/m/n/12 + \OML/cmm/m/it/12 c; n \OT1/cmr/m/n/12 : \O ML/cmm/m/it/12 b \OT1/cmr/m/n/12 + \OML/cmm/m/it/12 d; k \OT1/cmr/m/n/12 : \OML /cmm/m/it/12 a \OT1/cmr/m/n/12 + \OML/cmm/m/it/12 b; x\OT1/cmr/m/n/12 00 : \OML /cmm/m/it/12 a; lo \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 max\OT1/cmr/m/n/12 (0\OML /cmm/m/it/12 ; k \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 n\OT1/cmr/m/n/12 )\OML/c mm/m/it/12 ; hi \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 min\OT1/cmr/m/n/12 (\OML/cmm /m/it/12 k; m\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; support \OT1/cmr/m/n/12 : \OML /cmm/m/it/12 makelist\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 i; i; lo; hi\OT1/cmr/m/n /12 )\OML/cmm/m/it/12 ; logdc \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 create[]ist\OT 1/cmr/m/n/12 (\OML/cmm/m/it/12 log\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 pdf[]yperge ometric\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 i; m; n; k\OT1/cmr/m/n/12 ))\OML/cmm/m /it/12 ; i; support\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; doubleEps \OT1/cmr/m/n/1 2 :
Overfull \hbox (39.36833pt too wide) in paragraph at lines 143--239 \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 :\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 e \OMS/cms y/m/n/12 ^^@ \OT1/cmr/m/n/12 10\OML/cmm/m/it/12 ; dnhyper \OT1/cmr/m/n/12 : \OM L/cmm/m/it/12 lambda\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 ncp\OT1/cmr/m/n/12 ]\OML /cmm/m/it/12 ; block\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 maxd; sumd\OT1/cmr/m/n/1 2 2\OML/cmm/m/it/12 ; d; d\OT1/cmr/m/n/12 2]\OML/cmm/m/it/12 ; d \OT1/cmr/m/n/1 2 : \OML/cmm/m/it/12 makelist\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 logdc\OT1/cmr/m/ n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ] + \OML/cmm/m/it/12 log\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp\OT1/cmr/m/n/12 ) \OMS/cmsy/m/n/12 ^^C \OML/cmm/m/it/12 s upport\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; i; \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 ; length\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 lo gdc\OT1/cmr/m/n/12 ))\OML/cmm/m/it/12 ; maxd \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 apply\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 max; d\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; d\OT1/cmr/m/n/12 2 :
Overfull \hbox (265.6928pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 create[]ist\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 exp\OT1/cmr/m/n/1 2 (\OML/cmm/m/it/12 di \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 maxd\OT1/cmr/m/n/1 2 )\OML/cmm/m/it/12 ; di; d\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; sumd\OT1/cmr/m/n /12 2 : \OML/cmm/m/it/12 lsum\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 i; i; d\OT1/cmr/ m/n/12 2)\OML/cmm/m/it/12 ; create[]ist\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 d\OT1/ cmr/m/n/12 2\OML/cmm/m/it/12 i=sumd\OT1/cmr/m/n/12 2\OML/cmm/m/it/12 ; d\OT1/cm r/m/n/12 2\OML/cmm/m/it/12 i; d\OT1/cmr/m/n/12 2)))\OML/cmm/m/it/12 ; mnhyper \ OT1/cmr/m/n/12 : \OML/cmm/m/it/12 lambda\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 ncp\ OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; block\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 d; v alue\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; ifequal\OT1/cmr/m/n/12 (\OML/cmm/m/it/1 2 ncp; \OT1/cmr/m/n/12 0)\OML/cmm/m/it/12 thenloelseifequal\OT1/cmr/m/n/12 (\OM L/cmm/m/it/12 ncp; inf\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 thenhielseblock\OT1/cmr /m/n/12 ([]\OML/cmm/m/it/12 ; d \OT1/cmr/m/n/12 :
Overfull \hbox (64.61157pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 dnhyper\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp\OT1/cmr/m/n/12 )\ OML/cmm/m/it/12 ; sum\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 support\OT1/cmr/m/n/12 [ \OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ] \OMS/cmsy/m/n/12 ^^C \OML/cmm/m/it/12 d\OT1 /cmr/m/n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; i; \OT1/cmr/ m/n/12 1\OML/cmm/m/it/12 ; length\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 d\OT1/cmr/m/ n/12 )))))\OML/cmm/m/it/12 ; pnhyper \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 lambda\ OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 q; ncp; upperTail\OT1/cmr/m/n/12 ]\OML/cmm/m/ it/12 ; ifequal\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp; \OT1/cmr/m/n/12 1)\OML/cm m/m/it/12 then\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifupperTailthen\OT1/cmr/m/n/12 1 \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 cdf[]ypergeometric\OT1/cmr/m/n/12 (\OML /cmm/m/it/12 q \OMS/cmsy/m/n/12 ^^@
Overfull \hbox (236.55396pt too wide) in paragraph at lines 143--239 \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 ; m; n; k\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 el secdf[]ypergeometric\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 q; m; n; k\OT1/cmr/m/n/12 ))\OML/cmm/m/it/12 elseifequal\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp; \OT1/cmr/ m/n/12 0)\OML/cmm/m/it/12 then\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifupperTailthen \OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifq <\OT1/cmr/m/n/12 = \OML/cmm/m/it/12 lothe n\OT1/cmr/m/n/12 1\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 0)\OML/cmm/m/it/12 else\ OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifq >\OT1/cmr/m/n/12 = \OML/cmm/m/it/12 lothen \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 0))\OML/cmm/m/it/12 elsei fequal\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp; inf\OT1/cmr/m/n/12 )\OML/cmm/m/it/ 12 then\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifupperTailthen\OT1/cmr/m/n/12 (\OML/c mm/m/it/12 ifq <\OT1/cmr/m/n/12 =
Overfull \hbox (130.16159pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 hithen\OT1/cmr/m/n/12 1\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 0) \OML/cmm/m/it/12 else\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ifq >\OT1/cmr/m/n/12 = \ OML/cmm/m/it/12 hithen\OT1/cmr/m/n/12 1\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 0)) \OML/cmm/m/it/12 elseblock\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 d\OT1/cmr/m/n/12 ] \OML/cmm/m/it/12 ; d \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 dnhyper\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 ncp\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; ifupperTailthensum\OT1 /cmr/m/n/12 (\OML/cmm/m/it/12 ifsupport\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 i\OT1/ cmr/m/n/12 ] \OML/cmm/m/it/12 >\OT1/cmr/m/n/12 = \OML/cmm/m/it/12 qthend\OT1/cm r/m/n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 else\OT1/cmr/m/n/ 12 0\OML/cmm/m/it/12 ; i; \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 ; length\OT1/cmr/m/ n/12 (\OML/cmm/m/it/12 d\OT1/cmr/m/n/12 ))\OML/cmm/m/it/12 elsesum\OT1/cmr/m/n/ 12 (\OML/cmm/m/it/12 ifsupport\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/1 2 ] \OML/cmm/m/it/12 <\OT1/cmr/m/n/12 =
Overfull \hbox (181.71327pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 qthend\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ]\OML /cmm/m/it/12 else\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 ; i; \OT1/cmr/m/n/12 1\OML/c mm/m/it/12 ; length\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 d\OT1/cmr/m/n/12 ))))\OML/ cmm/m/it/12 ; ncpL \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 lambda\OT1/cmr/m/n/12 ([\ OML/cmm/m/it/12 alpha\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; ifequal\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; lo\OT1/cmr/m/n/12 )\OM L/cmm/m/it/12 then\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 elseblock\OT1/cmr/m/n/12 ([ \OML/cmm/m/it/12 p\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; p \OT1/cmr/m/n/12 : \OML/ cmm/m/it/12 pnhyper\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cm m/m/it/12 ; \OT1/cmr/m/n/12 1\OML/cmm/m/it/12 ; true\OT1/cmr/m/n/12 )\OML/cmm/m /it/12 ; ifp > alphathenfind[]oot\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 lambda\OT1/c mr/m/n/12 ([\OML/cmm/m/it/12 y\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; pnhyper\OT1/c mr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; y; true\OT1/c mr/m/n/12 ) \OMS/cmsy/m/n/12 ^^@
Overfull \hbox (3.40475pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 alpha\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; \OT1/cmr/m/n/12 0\OML /cmm/m/it/12 ; \OT1/cmr/m/n/12 1)\OML/cmm/m/it/12 elseifp < alphathen\OT1/cmr/m /n/12 1\OML/cmm/m/it/12 =find[]oot\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 lambda\OT1/ cmr/m/n/12 ([\OML/cmm/m/it/12 y\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; pnhyper\OT1/ cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; \OT1/cmr/m/n /12 1\OML/cmm/m/it/12 =y; true\OT1/cmr/m/n/12 ) \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m /it/12 alpha\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; doubleEps; \OT1/cmr/m/n/12 1)\O ML/cmm/m/it/12 else\OT1/cmr/m/n/12 1))\OML/cmm/m/it/12 ; ncpU \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 lambda\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 alpha\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; ifx\OT1/cmr/m/n/12 00 =
Overfull \hbox (274.29704pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 hitheninfelseblock\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 p\OT1/cmr /m/n/12 ]\OML/cmm/m/it/12 ; p \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 pnhyper\OT1/cm r/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; \OT1/cmr/m/n/1 2 1\OML/cmm/m/it/12 ; false\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; ifp < alphathenf ind[]oot\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 lambda\OT1/cmr/m/n/12 ([\OML/cmm/m/it /12 y\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; pnhyper\OT1/cmr/m/n/12 (\OML/cmm/m/it/ 12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; y; false\OT1/cmr/m/n/12 ) \OMS/cmsy/m/ n/12 ^^@ \OML/cmm/m/it/12 alpha\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; \OT1/cmr/m/n /12 0\OML/cmm/m/it/12 ; \OT1/cmr/m/n/12 1)\OML/cmm/m/it/12 elseifp > alphathen\ OT1/cmr/m/n/12 1\OML/cmm/m/it/12 =find[]oot\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 la mbda\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 y\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; pnh yper\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; \OT 1/cmr/m/n/12 1\OML/cmm/m/it/12 =y; false\OT1/cmr/m/n/12 ) \OMS/cmsy/m/n/12 ^^@
Overfull \hbox (190.69191pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 alpha\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; doubleEps; \OT1/cmr/m /n/12 1)\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 1))\OML/cmm/m/it/12 ; pvalue \OT1/ cmr/m/n/12 : (\OML/cmm/m/it/12 ifalternative \OT1/cmr/m/n/12 =[] \OML/cmm/m/it/ 12 lessthenpnhyper\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm /m/it/12 ; oddsratio; false\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 elseifalternative \OT1/cmr/m/n/12 =[] \OML/cmm/m/it/12 greaterthenpnhyper\OT1/cmr/m/n/12 (\OML/cm m/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; oddsratio; true\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 elseifalternative \OT1/cmr/m/n/12 =[]
Overfull \hbox (84.44579pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 twosidedthenblock\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 relerr; ds tar; dn\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 ; relErr \OT1/cmr/m/n/12 : 1 + 1\OML/c mm/m/it/12 :\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 e \OMS/cmsy/m/n/12 ^^@ \OT1/cmr/m /n/12 7\OML/cmm/m/it/12 ; dn \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 dnhyper\OT1/cmr /m/n/12 (\OML/cmm/m/it/12 oddsratio\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 ; dstar \O T1/cmr/m/n/12 : \OML/cmm/m/it/12 dn\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 x\OT1/cmr/ m/n/12 00 \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 lo \OT1/cmr/m/n/12 + 1] \OMS/cm sy/m/n/12 ^^C \OML/cmm/m/it/12 relErr; lsum\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 if di < dstarthendielse\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 ; di; dn\OT1/cmr/m/n/12 ) )\OML/cmm/m/it/12 else \OMS/cmsy/m/n/12 ^^@
Overfull \hbox (20.20848pt too wide) in paragraph at lines 143--239 \OT1/cmr/m/n/12 1)\OML/cmm/m/it/12 ; cInterval \OT1/cmr/m/n/12 : (\OML/cmm/m/it /12 ifalternative \OT1/cmr/m/n/12 =[] \OML/cmm/m/it/12 lessthen\OT1/cmr/m/n/12 [0\OML/cmm/m/it/12 ; ncpU\OT1/cmr/m/n/12 (1 \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/ 12 confLevel\OT1/cmr/m/n/12 )]\OML/cmm/m/it/12 elseifalternative \OT1/cmr/m/n/1 2 =[] \OML/cmm/m/it/12 greaterthen\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 ncpL\OT1/cm r/m/n/12 (1 \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 confLevel\OT1/cmr/m/n/12 )\OM L/cmm/m/it/12 ; inf\OT1/cmr/m/n/12 ]\OML/cmm/m/it/12 elseifalternative \OT1/cmr /m/n/12 =[]
Overfull \hbox (72.89813pt too wide) in paragraph at lines 143--239 \OML/cmm/m/it/12 twosidedthen\OT1/cmr/m/n/12 [\OML/cmm/m/it/12 ncpL\OT1/cmr/m/n /12 ((1 \OMS/cmsy/m/n/12 ^^@ \OML/cmm/m/it/12 confLevel\OT1/cmr/m/n/12 )\OML/cm m/m/it/12 =\OT1/cmr/m/n/12 2)\OML/cmm/m/it/12 ; ncpU\OT1/cmr/m/n/12 ((1 \OMS/cm sy/m/n/12 ^^@ \OML/cmm/m/it/12 confLevel\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 =\OT1 /cmr/m/n/12 2)]\OML/cmm/m/it/12 else\OT1/cmr/m/n/12 [\OMS/cmsy/m/n/12 ^^@\OT1/c mr/m/n/12 1\OML/cmm/m/it/12 ; \OMS/cmsy/m/n/12 ^^@\OT1/cmr/m/n/12 1])\OML/cmm/m /it/12 ; mle \OT1/cmr/m/n/12 : (\OML/cmm/m/it/12 ifequal\OT1/cmr/m/n/12 (\OML/c mm/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; lo\OT1/cmr/m/n/12 )\OML/cmm/m/ it/12 then\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 elseifequal\OT1/cmr/m/n/12 (\OML/cm m/m/it/12 x\OT1/cmr/m/n/12 00\OML/cmm/m/it/12 ; hi\OT1/cmr/m/n/12 )\OML/cmm/m/i t/12 theninfelseblock\OT1/cmr/m/n/12 ([\OML/cmm/m/it/12 mu\OT1/cmr/m/n/12 ]\OML /cmm/m/it/12 ; mu \OT1/cmr/m/n/12 : [3]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.240 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.242 \end{maxima} \newpage [4]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.243 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.245 \end{maxima} \newpage [5]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.246 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.248 \end{maxima} \newpage [6]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.249 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.253 \end{maxima} \newpage [7]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.254 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.263 \end{maxima} \newpage
Overfull \hbox (10.9121pt too wide) in paragraph at lines 255--263 []\T1/cmr/m/n/12 testSet() := block([fit1], fit1 : fish-ertest(10,10,10,20,'two sided,1,0.95), [test-Tol-er-ance(fit1[1], 0.2575, 1.0e-3), test-Tol-er-ance(fit 1[2][1],0.5383996, [8]
LaTeX Error: Environment maxima undefined.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.264 \begin{maxima}
LaTeX Error: \begin{document} ended by \end{maxima}.
See the LaTeX manual or LaTeX Companion for explanation. Type H <return> for immediate help. ...
l.266 \end{maxima} \newpage [9] [10] (./8783698267322362323-16.0px.aux) ) (see the transcript file for additional information) Output written on 8783698267322362323-16.0px.dvi (10 pages, 8592 bytes). Transcript written on 8783698267322362323-16.0px.log.