R/U-E1E1.C.R

Defines functions findmax p.ll.C p.est.E1E1.C p.estFUN.E1E1.C llsearch.E1E1.C

# Linearizable U-E1E1
# Common variance for both segments

llsearch.E1E1.C <- function(x, y, n, jlo, jhi)
{
	fj <- matrix(0, n)
	fxy <- matrix(0, jhi - jlo + 1)
	
	jgrid <- expand.grid(jlo:jhi)
	k.ll <- apply(jgrid, 1, p.estFUN.E1E1.C, x = x, y = y, n = n)
	
	fxy <- matrix(k.ll, nrow = jhi-jlo+1)
	rownames(fxy) <- jlo:jhi
  
	z <- findmax(fxy)
	jcrit <- z$imax + jlo - 1
	list(jhat = jcrit, value = max(fxy))
}

#  Function for deriving the ML estimates of the change-points problem.

p.estFUN.E1E1.C <- function(j, x, y, n){
  a <- p.est.E1E1.C(x,y,n,j)
  s2 <- a$sigma2
  return(p.ll.C(n, j, s2))
}

p.est.E1E1.C <- function(x,y,n,j){
 xa <- x[1:j]
 ya <- y[1:j]
 jp1 <- j+1
 xb <- x[jp1:n]
 yb <- y[jp1:n]
 g1 <- lm(ya ~ exp(xa))
 g2 <- lm(yb ~ exp(xb))
 beta <-c(g1$coef[1],g1$coef[2],g2$coef[1],g2$coef[2])
 s2 <- (sum((ya-g1$fit)^2)+sum((yb-g2$fit)^2))/n
list(a0=beta[1],a1=beta[2],b0=beta[3],b1=beta[4],sigma2=s2,xj=x[j])
}

#  Function to compute the log-likelihood of the change-point problem

p.ll.C <- function(n, j, s2){
 q1 <- n * log(sqrt(2 * pi))
 q2 <- 0.5 * n  * (1 + log(s2))
 - (q1 + q2)
}

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)
}

Try the vrcp package in your browser

Any scripts or data that you put into this service are public.

vrcp documentation built on May 29, 2017, 3:03 p.m.