#######Pearson's chisquare test for general 2 by J table x######
	mychisq=function(x){  
		n=sum(x)
		apply(x,1,sum)->row.sum
		apply(x,2,sum)->col.sum
		e=row.sum%*%t(col.sum)/n
		#e[e==0]=0.5
	sum((x-e)^2/e)
	}
##########

#####function to generate multinomial random number###
rmn=function(n,size,prob){
       U=matrix( runif(n*size),n, size)
       f.local=function(u, prob){
       k=length(prob)
       numb=1:k
       cp=c(0, cumsum(prob))
       for (i in 1:k){
               numb[i]=sum(cp[i]<u & u<=cp[i+1])
       }
       numb
        }
        apply(U, 1, f.local, prob=prob)
 }
######################################################



####Additive effect score function for K-allelic loci###
score<-function(K=3){
s=NULL
for(m in 1:K){
for(i in 1:K){
j=i
while(j<=K)
   {c(s,as.numeric(i==m)+as.numeric(j==m))->s
    j=j+1}
}
}
matrix(s,ncol=K)->score
score[,-K]}

###########################################################



## Cochran-Armitage's trend test: for K-allelic loci###
### score matrix x is preassigned as additive effect######

   CATK=function(data,K){
      x=score(K=K)
      N=sum(data)
     
	apply(data, 1, sum)-> rs
	apply(data,2, sum)->n
  	p.margin=n/N
	pr=data/rs
      phi=rs[1]/N
      sigma=N*(t(x)%*%diag(n)%*% x )-t(x)%*%n%*%t(n)%*%x
      diag(rep(1,K-1))->I
      solve(sigma,I)->sigmainv

      delta=t(x)%*%(data[1,]-phi*n)
      		
	test=N^3/(rs[1]*(N-rs[1]))*t(delta)%*%sigmainv%*%delta	
	test=as.numeric(test)
	test
      } #end of CATK with df=K-1 chisq with df=K(K+1)/2-1   


####generate contigency table with  J by K score matrix x#####

genK.data =function(n1=500, n0=500,alpha=-2,beta=rep(0,K-1),p,K,x=matrix(0,nrow=K*       (K+1)/2,ncol=K-1)){ 
     	## n1: numb case; n0:control;
        ## beta: (k-1)*1 vector
        ##p: the frequency of allele (M1,...MK)
    
        ### generate the genotype frequency in population###
        pgeno=NULL
        for(i in 1:K){
         j=i
         while(j<=K){
         if(i==j) pgeno=c(pgeno,p[i]^2)
         else     pgeno=c(pgeno,2*p[i]*p[j])
         j=j+1}  
        }
       ####generate penertrance for each genotype###

        f=1/(1+exp(-alpha-x%*%beta))
        p.case=f*pgeno
        p.cont=(1-f)*pgeno
        p.case=p.case/sum(p.case)
        p.cont=p.cont/sum(p.cont)
    	data = cbind(rmn(1, n1, p.case),rmn(1, n0, p.cont))
	data=t(data)
	dimnames(data)=list(c("case", "cont"), 1:(K*(K+1)/2))
	data


       } #end of genK.data function



###CATK2: Cochran-Armitage trend test with score vector (x1,x2,...xj)###

CATK2<-function(data,x){
apply(data, 1, sum)-> ns
apply(data,2, sum)->a
p.margin=a/sum(ns)
pr=data/ns
      delta=x%*%(pr[1,]-pr[2,])
sigma=(1/ns[1]+1/ns[2]) *( x^2%*%p.margin-(x%*%p.margin)^2  )   
test=delta^2/sigma
test=as.numeric(test)
return(test)
      } #end of CATK2 df=1



#####generate general 2*J table with score vector x: J=4###
genJ.data<-function(n1=n1,n2=n2,n3=n3,n4=n4,alpha=-2,beta=0,x=seq(from=1,to=J,by=1)){

f=1/(1+exp(-alpha- x*beta))
 p1=c(f[1],1-f[1])
 p2=c(f[2],1-f[2])
 p3=c(f[3],1-f[3])
 p4=c(f[4],1-f[4]) 

data=cbind(rmn(1,n1,p1),rmn(1,n2,p2),rmn(1,n3,p3),rmn(1,n4,p4))

dimnames(data)=list(c("Y=1", "Y=0"), 1:J)
data
}