R help - Feb 2005 - efficient R code

Last Friday, Gregory Chaitin (http://www.umcs.maine.edu/~chaitin/lm.html) 
mentioned that there can be no proof that a given code is the shortest for 
a problem, even within a language. Still, the script below, a replacement 
of the "TDT", one of the most frequently used tests in genetics 
(http://mustat.rockefeller.edu under "downloads") may get close. It 
contains a few additional bytes for clarity, as in (2^1) for 2, but, 
otherwise, I don't think one could make this much shorter, especially the 
part that does the exact distribution.

I'm looking forward to comments on improving the programming efficiency for 
this problem. (The "return(...)" seems to be necessary in R only.)

Knut

#--------------------------------------------------------------------------
# asymp.SMN.pvalue(pP,qP,pX,qX,pQ,qQ)
# exact.SMN.pvalue(pP,qP,pX,qX,pQ,qQ)
#--------------------------------------------------------------------------
# pP,qP = number of PP,PQ children of PP~PQ parents
# pX,qX = number of PP,QQ children of PQ~PQ parents
# pQ,qQ = number of PQ,QQ children of PQ~QQ parents
#--------------------------------------------------------------------------

O3 <- function(X1,X2,X3,Op) matrix(outer(outer(X1,X2,Op),X3,Op))
b0 <- function(n) if (n==0) 1 else dbinom(0:n, n, .5)

Dst <- function(nP,   nX,   nQ   )
O3(b0(nP),(2^0)*b0(nX),b0(nQ),"*")
Est <- function(pP,qP,pX,qX,pQ,qQ) O3(pP-qP,
(2^1)*(pX-qX),pQ-qQ,"+") # (1)
Var <- function(pP,qP,pX,qX,pQ,qQ) O3(pP+qP,
(2^2)*(pX+qX),pQ+qQ,"+") # (2)

asymp.SMN.pvalue <- function(...)  1-pchisq( Est(...)^2/Var(...) ,1)  # (3)

exact.SMN.pvalue <- function(pP,qP,pX,qX,pQ,qQ)
{
         tb <- cbind(
                 Dst(nP<-pP+qP, nX<-pX+qX, nQ<-pQ+qQ),
                 Est(0:nP,nP:0, 0:nX,nX:0, 0:nQ,nQ:0)^2)

         return(1-sum(tb[tb[,2]<c(Est(pP,qP,pX,qX,pQ,qQ)^2),1]))
}
#--------------------------------------------------------------------------
# Note: The numbers in parentheses in inline comments refer to the equation
#       numbers in
#       Wittkowski KM, Liu X
#       A Statistically Valid Alternative to the TDT
#       Hum Hered 2002, 54:157-164; 2004, 58:60-61 (discussion)
#--------------------------------------------------------------------------


Knut M. Wittkowski, PhD,DSc
------------------------------------------
The Rockefeller University, GCRC
Experimental Design and Biostatistics
1230 York Ave #121B, Box 322, NY,NY 10021
+1(212)327-7175, +1(212)327-8450 (Fax)
kmw at rockefeller.edu
http://www.rucares.org/clinicalresearch/dept/biometry/

On Wed, 9 Feb 2005, Knut M. Wittkowski wrote:
> Last Friday, Gregory Chaitin (http://www.umcs.maine.edu/~chaitin/lm.html) 
> mentioned that there can be no proof that a given code is the shortest for
a
> problem, even within a language.
I presume he said something more like `for some problems there is no 
proof' or `there is no procedure for constructing a proof'

You can prove that function(x) x+1 is the shortest function for computing 
x+1 in R by listing all the shorter functions and verifying that none of 
them computes x+1, both of which tasks can be done fairly fast.  This 
approach doesn't generalise because you can't give a general way to 
perform the second step: checking that none of the smaller programs works.

I would expect that the TDT falls in the category of problems checkable by 
enumeration.

>				Still, the script below, a replacement of 
> the "TDT", one of the most frequently used tests in genetics 
> (http://mustat.rockefeller.edu under "downloads") may get close.
It contains
> a few additional bytes for clarity, as in (2^1) for 2, but, otherwise, I 
> don't think one could make this much shorter, especially the part that
does
> the exact distribution.
>
> I'm looking forward to comments on improving the programming efficiency
for
> this problem. (The "return(...)" seems to be necessary in R
only.)
>
The return() should not be necessary in R.

 	-thomas

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle