swh:1:snp:ffdd0a7d2c8ea15ad41d45b3b178f668bd942287
Tip revision: 65dce028095b06ee5ed19673e3045548876c41d3 authored by Derek Young on 30 December 2007, 00:00:00 UTC
version 0.3.3
version 0.3.3
Tip revision: 65dce02
regmixchpt.R
regmixEM.chgpt <- function (y, x, lambda = NULL, gamma = NULL, beta = NULL, sigma = NULL, k = 2, T = 3, t=NULL,
epsilon = 1e-08, maxit = 10000, verb=FALSE)
{
x = cbind(1, x)
n <- length(y)
p <- ncol(x)
if(is.null(t)) {
stop(paste("You must specify a changepoint value!", "\n")) }
tmp<-regmix.chgpt.init(y=y, x=x, lambda=lambda, beta=beta, gamma=gamma, sigma=sigma, t=t, k=k)
lambda<-tmp$lambda
beta<-as.vector(tmp$beta)
gamma<-as.vector(tmp$gamma)
sigma<-tmp$sigma
diff <- 1
iter <- 0
xbeta <- x%*%beta
comp2 <- lambda[2]*dnorm(y,mean=xbeta,sd=sigma[2])
obsloglik <- NULL
seq <- T:(n-(T+1))
for(i in seq){
x.t <- (x[,2]-x[i,2])*(x[,2]>x[i,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
obsloglik[i-(T-1)] <- sum(log(apply(comp, 1, sum)))
}
if(sum(obsloglik==-Inf)==length(obsloglik)){
cutpt <- seq[floor(mean(c(T,(n-T-1))))]}
else cutpt <- seq[which.max(obsloglik)]
x.t <- (x[,2]-x[cutpt,2])*(x[,2]>x[cutpt,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
compsum <- apply(comp,1,sum)
obsloglik <- max(obsloglik,na.rm=T)
ll<-obsloglik
restarts<-0
while (diff > epsilon && iter < maxit) {
z <- comp/compsum
z[,k]=1-apply(as.matrix(z[,(1:(k-1))]),1,sum)
lambda.new <- apply(z,2,mean)
gamma.new <- try(solve(t(cbind(x,x.t)) %*% (as.vector(z[,1]) *
cbind(x,x.t))) %*% (t(cbind(x,x.t)) %*% (as.vector(z[,1]) * y)),silent=TRUE)
beta.new <- try(solve(t(x) %*% (as.vector(z[,2]) *
x)) %*% (t(x) %*% (as.vector(z[,2]) * y)),silent=TRUE)
if(sum(sigma < 1e-08)>0 || sum(lambda < 1e-08)>0 || class(gamma.new)=="try-error" || class(beta.new)=="try-error"
|| sum(is.na(gamma.new)>0) || sum(is.na(beta.new)>0)){
cat("Need new starting values due to singularity...", "\n")
restarts <- restarts + 1
if(restarts>15) stop("Too many tries!")
tmp<-regmix.chgpt.init(y=y, x=x, t=t, k=k)
lambda<-tmp$lambda
beta<-tmp$beta
gamma<-tmp$gamma
sigma<-tmp$sigma
diff <- 1
iter <- 0
xbeta <- x%*%beta
comp2 <- lambda[2]*dnorm(y,mean=xbeta,sd=sigma[2])
obsloglik <- NULL
seq <- T:(n-(T+1))
for(i in seq){
x.t <- (x[,2]-x[i,2])*(x[,2]>x[i,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
obsloglik[i-(T-1)] <- sum(log(apply(comp, 1, sum)))
}
if(sum(obsloglik==-Inf)==length(obsloglik)){
cutpt <- seq[floor(mean(c(T,(n-T-1))))]}
else cutpt <- seq[which.max(obsloglik)]
cutpt <- seq[which.max(obsloglik)]
x.t <- (x[,2]-x[cutpt,2])*(x[,2]>x[cutpt,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
compsum <- apply(comp,1,sum)
obsloglik <- max(obsloglik,na.rm=T)
ll<-obsloglik
} else{
sigma.new <- NULL
sigma.new[1] <- sqrt(sum(as.vector(z[,1])*(y-cbind(x,x.t)%*%gamma.new)^2)/sum(z[,1]))
sigma.new[2] <- sqrt(sum(as.vector(z[,2])*(y-x%*%beta.new)^2)/sum(z[,2]))
gamma <- gamma.new
beta <- beta.new
sigma <- sigma.new
lambda <- lambda.new
xbeta <- x%*%beta
comp2 <- lambda[2]*dnorm(y,mean=xbeta,sd=sigma[2])
newobsloglik <- NULL
for(i in seq){
x.t <- (x[,2]-x[i,2])*(x[,2]>x[i,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
newobsloglik[i-(T-1)] <- sum(log(apply(comp, 1, sum)))
}
if(sum(newobsloglik==-Inf)==length(newobsloglik)){
cutpt <- seq[floor(mean(c(T,(n-T-1))))]}
else cutpt <- seq[which.max(newobsloglik)]
x.t <- (x[,2]-x[cutpt,2])*(x[,2]>x[cutpt,2])
xgamma <- cbind(x,x.t)%*%t(t(gamma))
comp1 <- lambda[1]*dnorm(y,mean=xgamma,sd=sigma[1])
comp <- cbind(comp1,comp2)
compsum <- apply(comp,1,sum)
newobsloglik <- max(newobsloglik,na.rm=T)
diff <- newobsloglik-obsloglik
obsloglik <- newobsloglik
ll <- c(ll, obsloglik)
iter<-iter+1
if (verb) {
cat("iteration=", iter, "diff=", diff, "log-likelihood",
obsloglik, "\n")
}
}
}
if (iter == maxit) {
cat("WARNING! NOT CONVERGENT!", "\n")
}
cat("number of iterations=", iter, "\n")
a=list(lambda = lambda, gamma = gamma, beta = beta, sigma = sigma, cutpoint=cutpt, loglik = obsloglik,
posterior = z, all.loglik=ll, restarts=restarts, ft="regmixEM.chgpt")
class(a) = "mixEM"
a
}