# source("http://stat-help.net/Goodman-Biometrika-61-1974.txt")
source("http://stat-help.net/model-H1.txt")
# # Goodman, Biometrika(1974) 61, 2, p.215 Table 1
# Call function to generate the model H_1
# and plot some output

starts=4; iters=120; # set.seed(10)
# set.seed(11)
res=get.g1(starts, iters)
# jpeg("Goodman-first-model-2-latent-classes.jpeg")
X11()
par(mar=c(5,5,6,2), oma=c(2,2,2,2))
plot(res$bigmat[1,], xlim=c(0, ncol(res$bigmat)), ylim=c(0,1), 
	ylab=expression(hat(pi)[1]^X),
	main=expression(paste(hat(pi)[1]^X," coverges to 0.279")),
	cex.main=1.5, cex.lab=1.3, pch=16, cex=0.4)
for(i in 2:starts) points(res$bigmat[9*(i-1)+1,], pch=16, col=i, cex=0.4)
mtext("Prob that an individual is in latent class 1",
	line=-0.5,outer=T, cex=1.5)

abline(h=0.279, lty=3, lwd=3, col="red")
axis(2, 0.279, " ", col="red")
mtext("0.279", 2, line=0.6, las=2, at=0.279, col="red")

# dev.off()


# Call function to generate the model H_4^{''}
# and plot some output
source("http://stat-help.net/model-H4-db-ABC-only.txt"); 
source("http://stat-help.net/model-H4-db-BCD-only.txt");
source("http://stat-help.net/model-H4-db-combine-estimates.txt")
source("http://stat-help.net/model-H4-db-final.txt")
# set.seed(21); 
#########################################################
# This first part generates deliberate starting values and
#  uses them to find final values
max.ii=3; starts=8; iters=100;
# Do this entire process a number (max.ii) of times
for(ii in 1:max.ii){
# estimate parameters in 2-class models for tables {ABC} and {BCD}
# completely restart with starts randomly generated sets of starting
# values and assume convergence after iters iterations 
res=get.g4.dp(starts, iters)
res2=get.g4.dp.b(starts, iters)

tmp=which(res$lik.rat.chsq==min(res$lik.rat.chsq, na.rm=T)) # 
tmp2=which(res2$lik.rat.chsq==min(res2$lik.rat.chsq, na.rm=T)) # 

tmp=res$bigmat[(7*(tmp-1)+1):(7*tmp), iters+1]
tmp2=res2$bigmat[(7*(tmp2-1)+1):(7*tmp2), iters+1]
# At this point, one has to take the most recently obtained estimates 
# and, using the hints given on p. 227 "Restriction (39) 
# in the 3-way table...",
# determine which class is the first class
coef.mat=matrix(c(1,1,0, 1,0,0, 0,1,1), nrow=3, ncol=3, byrow=T)
# A and C
t.A.and.C=solve(coef.mat, c(tmp[1], tmp2[1], 1-tmp2[1]))
t.A.and.C=sum(t.A.and.C>1.1 | t.A.and.C<(-0.1))==0
# A and D
t.A.and.D=solve(coef.mat, c(tmp[1], 1-tmp2[1], tmp2[1]))
t.A.and.D=sum(t.A.and.D>1.1 | t.A.and.D<(-0.1))==0
# B and C
t.B.and.C=solve(coef.mat, c(1-tmp[1], tmp2[1], 1-tmp2[1]))
t.B.and.C=sum(t.B.and.C>1.1 | t.B.and.C<(-0.1))==0
# B and D
t.B.and.D=solve(coef.mat, c(1-tmp[1], 1-tmp2[1], tmp2[1]))
t.B.and.D=sum(t.B.and.D>1.1 | t.B.and.D<(-0.1))==0

t.A.and.C.lik=Inf; t.A.and.D.lik=Inf; t.B.and.C.lik=Inf; t.B.and.D.lik=Inf
if(t.A.and.C) 
	t.A.and.C.lik=get.g4.dp.combine(c(tmp2[1], tmp[1]-tmp2[1], 1-tmp[1],
		tmp[c(2,2,3,4,4,5)], tmp2[c(4,5,5,6,7,7)]))  

if(t.A.and.D) 
	t.A.and.D.lik=get.g4.dp.combine(c(1-tmp2[1], tmp2[1]+tmp[1]-1, 1-tmp[1],
		tmp[c(2,2,3,4,4,5)], tmp2[c(5,4,4,7,6,6)]))  

if(t.B.and.C) 
	t.B.and.C.lik=get.g4.dp.combine(c(tmp2[1], 1-(tmp2[1]+tmp[1]), tmp[1],
		tmp[c(3,3,2,5,5,4)], tmp2[c(4,5,5,6,7,7)]))  

if(t.B.and.D) 
	t.B.and.D.lik=get.g4.dp.combine(c(1-tmp2[1], tmp2[1]-tmp[1], tmp[1],
		tmp[c(3,3,2,5,5,4)], tmp2[c(5,4,4,7,6,6)])) 

tvec=c(t.A.and.C.lik, t.A.and.D.lik, t.B.and.C.lik, t.B.and.D.lik)
tt=which(tvec==min(tvec[is.finite(tvec)]))
if(tt==1) start.params=c(tmp2[1], tmp[1]-tmp2[1], 1-tmp[1],
		tmp[c(2,2,3)],(tmp[c(4,4)]+tmp2[c(2,2)])/2,tmp[5], 
		tmp2[4],(tmp2[c(5,5)]+tmp[c(7,7)])/3,tmp2[c(6,7,7)]) 

if(tt==2) start.params=c(1-tmp2[1], tmp2[1]+tmp[1]-1, 1-tmp[1],
		tmp[c(2,2,3)],(tmp[c(4,4)]+tmp2[c(3,3)])/2, tmp[5], 
		tmp2[5], (tmp2[c(4,4)]+tmp[c(7,7)])/2,
		tmp2[c(7,6,6)])

if(tt==3) start.params=c(tmp2[1], 1-(tmp2[1]+tmp[1]), tmp[1],
		tmp[c(3,3,2)],(tmp[c(5,5)]+tmp2[c(2,2)])/2,tmp[4], 
		tmp2[4],(tmp[c(6,6)]+tmp2[c(5,5)])/2,
		tmp2[c(6,7,7)])

if(tt==4) start.params=c(1-tmp2[1], tmp2[1]-tmp[1], tmp[1],
		tmp[c(3,3,2)],(tmp[c(5,5)]+tmp2[c(3,3)])/2,tmp[4], 
		tmp2[5],(tmp2[c(4,4)]+tmp[c(6,6)])/2,
		tmp2[c(7,6,6)])
#################################################
# Each of the max.ii times that we pass through here,
# use deliberate values to find the best value
f.g4=get.g4.dp.final(start.params, iters=iters)

if(ii==1) {
# jpeg("Goodman-fourth-model-3-latent-classes-good-start-vals-80-starts.jpeg",
# width=700)
X11(width = 12)
par(mar=c(5,5,6,2), oma=c(2,2,2,2), mfcol=c(1,2))
	plot(f.g4$params[1,],ylab=expression(hat(pi)[1]^X),
		main="", ylim=c(0,1),
		cex.main=1.5, pch=16, cex=0.4, cex.lab=1.3)
	abline(h=0.2567286, v=iters/2)
	mtext("Prob that an individual is in latent class 1",
		line=0, outer=T, cex=1.5)
	mtext(expression(paste("Model ",H[4],"'', deliberate starting values")),
		line=-2.5, outer=T, cex=1.5, adj=0)
	}else{
	points(f.g4$params[1,], col=ii, pch=ii-1, cex=0.4)
	}
if(ii==max.ii) {
	abline(h=c(0.2567269, 0.712523), v=c(iters/4, iters/2))
	text(c(3*iters/4, 3*iters/4), c(0.2567269, 0.712523), pos=1, 
		c("y=0.257", "y=0.713"))
	}

} # end "for ii in 1:max.ii..."
####################################################################
# This next part finds final values starting from random, not deliberate,
# values. Do a total of random.starts random starts
# set.seed(200); 
random.starts=5
f.g4=get.g4.dp.final(runif(15), iters=iters)

plot(f.g4$params[1,], ylim=c(0,1),ylab=expression(hat(pi)[1]^X),
	main=" ",
	cex.main=1.5, cex.lab=1.3, pch=16, cex=0.4)
mtext(expression(paste("Model ",H[4],"'', random starts")),
	line=-2.5,outer=T, cex=1.5, adj=1)
abline(v=iters/2)

for(i in 1:(random.starts-1)){
	f.g4=get.g4.dp.final(runif(15), iters=iters)
	points(f.g4$params[1,], col=i+1, pch=i, cex=0.4)

}
abline(h=c(0.2567269, 0.712523), v=c(iters/4, iters/2))
text(c(3*iters/4, 3*iters/4), c(0.2567269, 0.712523), 
	pos=1, c("y=0.257", "y=0.713"))

# dev.off()

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