R/MethUnc.D.R

In hcp: Change Point Estimation for Regression with Heteroscedastic Data

# Unconstraint Method
# Different variances for all three segments


llsearch.D <-function(x, y, n, jlo, jhi, klo, khi,plot)
{            
	fjk <- matrix(0, n, n)
	fxy <- matrix(0, (jhi - jlo + 1), (khi - klo + 1))
## Yulei's edit to avoid using for-loop
	jkgrid <- expand.grid(jlo:jhi, klo:khi)
	res <- data.frame(j = jkgrid[,1],
	                  k = jkgrid[,2],
	                  k.ll = apply(jkgrid, 1, p.estFUN.D, x = x, 
	                               y = y, n = n))
	
	fxy <- matrix(res$k.ll, nrow = jhi-jlo+1, ncol = khi-klo+1)
	rownames(fxy) <- jlo:jhi
	colnames(fxy) <- klo:khi
if (plot == "TRUE") {
  jx<-jlo:jhi
  ky<-klo:khi
  persp(jx, ky, fxy, xlab = "j", ylab = "k", zlab = "LL(x,y,j,k)")
  title("Log-likelihood Surface")
}
	z <- findmax(fxy)
	jcrit <- z$imax + jlo - 1
	kcrit <- z$jmax + klo - 1
	list(jhat = jcrit, khat = kcrit, value = max(fxy))
}

p.estFUN.D <- function(jk, x, y, n){
  j = jk[1]
  k = jk[2]
  a <- p.est.D(x,y,n,j,k)
  s2 <- a$sigma2
  t2 <- a$tau2
  u2 <- a$u2
  return(p.ll.D(n, j, k, s2, t2, u2))
}

p.est.D <-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)
 beta <- c(g1$coef[1],g1$coef[2],g2$coef[1],g2$coef[2],g3$coef[1],g3$coef[2])
 s2 <- sum((ya-g1$fit)^2)/j  
 t2 <- sum((yb-g2$fit)^2)/(k-j)
## Yulei's edit: u2 is psi estimation for the third segment
 u2 <- sum((yc-g3$fit)^2)/(n-k)
 list(a0=beta[1],a1=beta[2],b0=beta[3],b1=beta[4],c0=beta[5],c1=beta[6],sigma2=s2,tau2=t2,u2=u2,xj=x[j],xk=x[k])
}

p.ll.D<-function(n, j, k, s2, t2, u2){
 q1 <- n * log(sqrt(2 * pi))
 q2 <- 0.5 * (j) * (1 + log(s2))
 q3 <- 0.5 * (k - j) * (1 + log(t2))
## Yulei's edit
q4 <- 0.5*(n-k)*(1+log(u2))
 - (q1 + q2 + q3 + q4)
}

## Yulei's edit to avoid using for-loop
findmax <-function(a)
  {
    maxa<-max(a)
    imax<- which(a==max(a),arr.ind=TRUE)[1]
    jmax<-which(a==max(a),arr.ind=TRUE)[2]
	list(imax = imax, jmax = jmax, value = maxa)
}