# source("G4-double-prime-final.pl")
get.g4.dp.final=function(start.params, iters=150, n=216){
# Goodman, Biometrika(1974) 61, 2, p.215 Table 1
# This function generates the model H_1
# 
 
 # initialize big results matrix and lik.rat.chsq
 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))
pi.ABCD.X=array(NA, dim=c(rep(2,4),3)) # This is conditional prob 
# that X=(certain value) given that A,B,C,D equal certain values

for(jj in 1:iters){
	
	# 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])
	
	# Eqn 11
	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.ABCD.X[i,j,k,l,t]=
				pi.ABCDX[i,j,k,l,t]/pi.ABCD[i,j,k,l]
	}
	pi.ABCD.X[i,j,k,l,1]=1-(sum(pi.ABCD.X[i,j,k,l,])-pi.ABCD.X[i,j,k,l,1])
	}}}}


	
	# Eqn 12
	pi.X[1]=sum(pi.ABCD.X[,,,,1]*p)
	pi.X[2]=sum(pi.ABCD.X[,,,,2]*p)
	pi.X[3]=1-sum(pi.X[-3])


	# Instead of Eqns 13a-13d, impose restrictions
	t=1
	pi.A.X[1,1:2]=sum(p[1,,,]*
		(pi.ABCD.X[1,,,,t]/pi.X[t]+pi.ABCD.X[1,,,,t+1]/pi.X[t+1]))/2

	pi.B.X[1,1:2]=sum(p[,1,,]*
		(pi.ABCD.X[,1,,,t]/pi.X[t]+pi.ABCD.X[,1,,,t+1]/pi.X[t+1]))/2

	
	pi.C.X[1,t]=sum(p[,,1,]*pi.ABCD.X[,,1,,t])/pi.X[t]
	pi.D.X[1,t]=sum(p[,,,1]*pi.ABCD.X[,,,1,t])/pi.X[t]
	
	t=3
	pi.A.X[1,3]=sum(p[1,,,]*pi.ABCD.X[1,,,,t])/pi.X[t]
	pi.B.X[1,3]=sum(p[,1,,]*pi.ABCD.X[,1,,,t])/pi.X[t]
	pi.C.X[1,2:3]=sum(p[,,1,]*
		(pi.ABCD.X[,,1,,t-1]/pi.X[t-1]+pi.ABCD.X[,,1,,t]/pi.X[t]))/2
	
	pi.D.X[1,2:3]=sum(p[,,,1]*
		(pi.ABCD.X[,,,1,t-1]/pi.X[t-1]+pi.ABCD.X[,,,1,t]/pi.X[t]))/2

	pi.A.X[2,]=1-pi.A.X[1,]
	pi.B.X[2,]=1-pi.B.X[1,]
	pi.C.X[2,]=1-pi.C.X[1,]
	pi.D.X[2,]=1-pi.D.X[1,]
	
		
# cat("tmp",tmp,"pi.X[2]",pi.X[2],"\n")
params=cbind(params, c(pi.X[1:2], pi.A.X[1,], pi.B.X[1,], 
		pi.C.X[1,], pi.D.X[1,]))
lik.rat.chsq=c(lik.rat.chsq, 2*sum(n*p*log(p/pi.ABCD)))
  } # Closes "for(jj in 1:iters){"

  
return(list(params=params, lik.rat.chsq=lik.rat.chsq, p=p))
} # Closes f'n get.g1