#  MethMMP.s  
#
#  S-Plus functions used in the
#  Maximization-Maximization-Posterior method

#  mmpiter.s
mmpiter <- function(x,y,n,jlo,jhi,klo,khi,sig0,tau0,rr,ss,vv,ww)
{
xi <- matrix(0,10,4)
xi[1,3] <- sig0
xi[1,4] <- tau0
for (iter in 2:10){
 a <- mmp(x,y,n,jlo,jhi,klo,khi,xi[iter-1,3],xi[iter-1,4],rr,ss,vv,ww)
 xi[iter,1] <- a$jnew
 xi[iter,2] <- a$knew
 xi[iter,3] <- a$signew
 xi[iter,4] <- a$taunew
 diff1 <- xi[iter,1]-xi[iter-1,1]
 diff2 <- xi[iter,2]-xi[iter-1,2]
 if (abs(diff1) < .1 && abs(diff2) < .1) break
 }
 list(jhat=a$jnew,khat=a$knew,sigma2hat=a$signew,tau2hat=a$taunew)
}

#  mmp.s
mmp <- function(x, y, n, jlo, jhi, klo, khi, sig, tau, rr, ss, vv, ww)
{
	# 1.  Given sigma^2 & tau^2, find j & k
	step1 <- newjk(x, y, n, jlo, jhi, klo, khi, sig, tau)
	jhat <- step1$jhat
	khat <- step1$khat

	# 2.  Given j & k, find sigma^2 & tau^2
	step2 <- newst(x, y, n, jhat, khat, rr, ss, vv, ww)
	sighat <- step2$sig2hat
	tauhat <- step2$tau2hat
	list(jnew = jhat, knew = khat, signew = sighat, taunew = tauhat)
}

#  newjk.s:  Maximize log[pi(j,k|...)]
newjk <- function(x, y, n, jlo, jhi, klo, khi, sigma2, tau2)
{
	c <- (-10^120)
	fjk <- matrix(c, n, n)
	
	for(j in jlo:jhi) {
		for(k in klo:khi) {
			fjk[j, k] <- lpjk(x, y, n, j, k, sigma2, tau2)
		}
	}
	z <- findmax(fjk)
	jcrit <- z$imax
	kcrit <- z$jmax
	list(jhat = jcrit, khat = kcrit, value = max(fjk))
}

#  newst.s:  Maximize pi(sigma^2,tau^2|...)
newst <- function(x, y, n, j, k, rr, ss, vv, ww)
{
	g <- rss(x, y, n, j, k)
	rss1 <- g$S1
	rss2 <- g$S2
	rss3 <- g$S3
	s2 <- (2 * ss + rss1 + rss3)/(n - k + j + 2 * rr + 2)
	t2 <- (2 * ww + rss2)/(k - j + 2 * vv + 2)
	list(sig2hat = s2, tau2hat = t2)
}

#  log p_jk
lpjk <- function(x, y, n, j, k, sigma2, tau2)
{
	g <- rss(x, y, n, j, k)
	rss1 <- g$S1
	rss2 <- g$S2
	rss3 <- g$S3
	x1 <- ((k - j)/2) * log(sigma2/tau2)
	x2 <- (rss1 + rss3)/(2 * sigma2)
	x3 <- rss2/(2 * tau2)
	x1 - x2 - x3
}

#  rss - Residual Sum of Squares
rss <- function(x, y, n, j, k)
{
	xa <- x[1:j]
	ya <- y[1:j]
	jp1 <- j + 1
	xb <- x[jp1:k]
	yb <- y[jp1:k]
	kp1 <- k + 1
	xc <- x[kp1:n]
	yc <- y[kp1:n]
	g1 <- lm(ya ~ xa)
	g2 <- lm(yb ~ xb)
	g3 <- lm(yc ~ xc)
	a0 <- g1$coef[1]
	a1 <- g1$coef[2]
	b0 <- g2$coef[1]
	b1 <- g2$coef[2]
	c0 <- g3$coef[1]
	c1 <- g3$coef[2]
	s1 <- sum((g1$res)^2)
	s2 <- sum((g2$res)^2)
	s3 <- sum((g3$res)^2)
	list(a0 = a0, a1 = a1, b0 = b0, b1 = b1, c0 = c0, c1 = c1, S1 = s1, S2 = s2, 
		S3 = s3)
}
findmax <- function(a)
{
	I <- dim(a)[1]
	J <- dim(a)[2]
	maxa <- max(a)
	for(i in 1:I) {
		for(j in 1:J) {
			temp <- a[i, j]
			if(temp == maxa)
				imax <- i
			if(temp == maxa)
				jmax <- j
		}
	}
	list(imax = imax, jmax = jmax, value = maxa)
}