# source("G4-double-prime-b.pl")
get.g4.dp.b=function(starts=5, iters=150){

  
 # initialize big results matrix and lik.rat.chsq
 bigmat=NULL; 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
 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,3))
# cat("OK3\n")
 for(j in 1:2){for(k in 1:2){for(l in 1:2){
	p[j,k,l]=
		sum(tmp[des[,2]==j &  des[,3]==k &  des[,4]==l ])
	
 }}}
rm(tmp) # Keep it clean fellas!
# cat("OK2\n")

for(ii in 1:starts){ 
	# This loop doesn't close 'til almost 
	# the end of the f'n
 # we need 14 independent params
 tmp=round(runif(7), 2)

# pi.X.0=tmp[1]; pi.X.1=1-tmp[1];
pi.X=c(tmp[1], 1-tmp[1])

# 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.B.X=rbind(tmp2<-c(tmp[2], tmp[3]), 
	1-tmp2)
pi.C.X=rbind(tmp2<-c(tmp[4], tmp[5]), 
	1-tmp2)
pi.D.X=rbind(tmp2<-c(tmp[6], tmp[7]), 
	1-tmp2)

params=tmp
# cat("OK1\n")

pi.BCDX=array(0, dim=rep(2,4))
pi.BCD=array(0, dim=rep(2,3))
pi.BCD.X=array(0, dim=rep(2,4)) # Goodman denotes this pi.ABCD.X

for(jj in 1:iters){
	# Eqn 10
	for(j in 1:2){for(k in 1:2){for(l in 1:2){
		for(t in 1:2){
			pi.BCDX[j,k,l,t]=pi.X[t]*pi.B.X[j,t]*
				pi.C.X[k,t]*pi.D.X[l,t]
	}}}}
	     	       
	# Eqn 9
	for(j in 1:2){for(k in 1:2){for(l in 1:2){
		pi.BCD[j,k,l]=sum(pi.BCDX[j,k,l,])
	}}}

	# Eqn 11
	for(j in 1:2){for(k in 1:2){for(l in 1:2){
		for(t in 1:2){
			pi.BCD.X[j,k,l,t]=
				pi.BCDX[j,k,l,t]/pi.BCD[j,k,l]
	}}}}

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


	# Eqns 13a-13d
	pi.B.X[1,1]=sum(p[1,,]*pi.BCD.X[1,,,1])/pi.X[1]
	pi.C.X[1,1]=sum(p[,1,]*pi.BCD.X[,1,,1])/pi.X[1]
	pi.D.X[1,1]=sum(p[,,1]*pi.BCD.X[,,1,1])/pi.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], pi.B.X[1,], pi.C.X[1,], 
		pi.D.X[1,]))
  } # Closes "for(jj in 1:150){"

bigmat=rbind(bigmat, params)
lik.rat.chsq=c(lik.rat.chsq, 2*sum(n*p*log(p/pi.BCD)))
} # Closes "for(ii in 1:starts){"
return(list(bigmat=bigmat, lik.rat.chsq=lik.rat.chsq, p=p))
} # Closes f'n get.g1