# source("G1c.pl")
get.g1=function(starts=5, iters=150){
# Goodman, Biometrika(1974) 61, 2, p.215 Table 1
# This function generates the model H_1
# 
 
 # 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 # 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!

for(ii in 1:starts){ 
	# This loop doesn't close 'til almost 
	# the end of the f'n
 # we need 9 params
 tmp=round(runif(9), 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.A.X=cbind(c(tmp[2], 1-tmp[2]), c(tmp[3], 1-tmp[3]))
pi.B.X=cbind(c(tmp[4], 1-tmp[4]), c(tmp[5], 1-tmp[5]))
pi.C.X=cbind(c(tmp[6], 1-tmp[6]), c(tmp[7], 1-tmp[7]))
pi.D.X=cbind(c(tmp[8], 1-tmp[8]), c(tmp[9], 1-tmp[9]))

params=tmp
# params=cbind(params, c(pi.X[1], pi.A.X[,1], pi.B.X[,1], 
# pi.C.X[,1], pi.D.X[,1]) 
rm(tmp)

pi.ABCDX=array(runif(32), dim=rep(2,5))
pi.ABCD=array(runif(16), dim=rep(2,4))
pi.ABCD.X=array(runif(32), dim=rep(2,5)) # This is conditional prob 
# that X=(certain value) given that A,B,C,D equal certain values

for(jj in 1:iters){
	# Eqn 12
	pi.X[1]=sum(pi.ABCD.X[,,,,1]*p)
	pi.X[2]=1-pi.X[1]


	# Eqns 13a-13d
	for(t in 1:2){
		pi.A.X[1,t]=sum(p[1,,,]*pi.ABCD.X[1,,,,t])/pi.X[t]
		pi.B.X[1,t]=sum(p[,1,,]*pi.ABCD.X[,1,,,t])/pi.X[t]
		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]
	}
	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,]
	
	# 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:2){
			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]
	}}}}}
	     	       
	# 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,])
	}}}}
	# But that is one more calculated than we need
	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:2){
			pi.ABCD.X[i,j,k,l,t]=
				pi.ABCDX[i,j,k,l,t]/pi.ABCD[i,j,k,l]
	}}}}}

	
# cat("tmp",tmp,"pi.X[2]",pi.X[2],"\n")
params=cbind(params, c(pi.X[1], pi.A.X[1,], pi.B.X[1,], 
		pi.C.X[1,], pi.D.X[1,]))
  } # Closes "for(jj in 1:iters){"

  
  # Now we have (no. params x no. iters) matrix of params
  # Let's allocate classes correctly
  # Goodman chooses (Tables 2 and 3) X=1 to be the smallest class
if(is.finite(pi.X[1])){ if(!(pi.X[1]<=pi.X[2])){ 
	params[1,]=1-params[1,]
	params[2:9,]=params[c(3,2,5,4,7,6,9,8),]
}} # close "if(is.finite(pi.X))..."
	# if this condition isn't true,
	# switch all the entries

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