# source("model-H4-db-combine-estimates.txt")
# source("G4-double-prime-combine-estimates.pl")
get.g4.dp.combine=function(start.params){

 lik.rat.chsq=NULL
 
 # make the m-way contingency table of observed data
 # Goodman, Biometrika(1974) 61, 2, p.215 Table 1
 des=1-expand.grid(0:1, 0:1, 0:1, 0:1)[,4:1]
 des=2-des # Note that, under Goodman's notation, 
 # p_{1,1,1,1}=P(A=1, B=1, C=1, D=1)=
 # P(A=+, B=+, C=+, D=+)=42/216=0.1944,
 # i.e. his negative "-" level corresponds to "2"
 tmp=c(42, 23, 6, 25, 6, 24, 7, 38, 1, 4, 1, 6, 2, 9, 2, 20)
 # n is number of subjects
 n=sum(tmp)
 des=cbind(des, tmp)
 tmp=tmp/sum(tmp)
 # initialize the observed matrix
 p=array(0, rep(2,4))

 for(i in 1:2){for(j in 1:2){for(k in 1:2){for(l in 1:2){
	p[i,j,k,l]=tmp[des[,1]==i &  des[,2]==j &  des[,3]==k &  des[,4]==l]
	
 }}}}
rm(tmp) # Keep it clean fellas!

 # we need 9 params
 tmp=start.params
# pi.X.0=tmp[1]; pi.X.1=1-tmp[1];
pi.X=c(tmp[1], tmp[2], tmp[3])

# pi.letter.X is the array prob A=1 given X=1, prob A=2 given X=1,
# prob A=1 given X=2, prob A=2 given X=2, 
# prob B=1 given X=1, prob B=2 given X=1, etc.
# first index is value of A, B, C or D
# second index is value of X
pi.A.X=cbind(c(tmp[4], 1-tmp[4]), c(tmp[5], 1-tmp[5]), c(tmp[6], 1-tmp[6]))
pi.B.X=cbind(c(tmp[7], 1-tmp[7]), c(tmp[8], 1-tmp[8]), c(tmp[9], 1-tmp[9]))
pi.C.X=cbind(c(tmp[10], 1-tmp[10]), c(tmp[11], 1-tmp[11]), c(tmp[12], 1-tmp[12]))
pi.D.X=cbind(c(tmp[13], 1-tmp[13]), c(tmp[14], 1-tmp[14]), c(tmp[15], 1-tmp[15]))

params=NULL
rm(tmp)

pi.ABCDX=array(NA, dim=c(rep(2,4),3))
pi.ABCD=array(NA, dim=rep(2,4))
	
	# Eqn 10
	for(i in 1:2){for(j in 1:2){for(k in 1:2){
		for(l in 1:2){for(t in 1:3){
			pi.ABCDX[i,j,k,l,t]=pi.X[t]*pi.A.X[i,t]*
				pi.B.X[j,t]*pi.C.X[k,t]*pi.D.X[l,t]
	}}}}}
	pi.ABCDX[1,1,1,1,1]=1-(sum(pi.ABCDX)-pi.ABCDX[1,1,1,1,1])
	       
	# Eqn 9
	for(i in 1:2){for(j in 1:2){for(k in 1:2){
		for(l in 1:2){
			pi.ABCD[i,j,k,l]=sum(pi.ABCDX[i,j,k,l,])
	}}}}
	pi.ABCD[1,1,1,1]=1-(sum(pi.ABCD)-pi.ABCD[1,1,1,1])
	

lik.rat.chsq=2*sum(n*p*log(p/pi.ABCD))

  
return(lik.rat.chsq)
} # Closes f'n get.g1