# BCPA/All BCPA Functions

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
GetRho <- function (x, t)
{
getL <- function(rho) {
dt <- diff(t)
s <- sd(x)
mu <- mean(x)
n <- length(x)
x.plus <- x[-1]
x.minus <- x[-length(x)]
Likelihood <- dnorm(x.plus, mean = mu + (rho^dt) * (x.minus -
mu), sd = s * sqrt(1 - rho^(2 * dt)))
logL <- sum(log(Likelihood))
if (!is.finite(logL))
logL <- -10^10
return(-logL)
}
o <- optimize(getL, lower = 0, upper = 1, tol = 1e-04)
return(c(o\$minimum, o\$objective))
}

GetLL <- function (x, t, mu, s, rho)
{
dt <- diff(t)
n <- length(x)
x.plus <- x[-1]
x.minus <- x[-length(x)]
Likelihood <- dnorm(x.plus, mean = mu + (rho^dt) * (x.minus -
mu), sd = s * sqrt(1 - rho^(2 * dt)))
LL <- -sum(log(Likelihood))
return(LL)
}

GetDoubleL <- function(x,t,tbreak)
{
x1 <- x[1:tbreak]
x2 <- x[tbreak:length(x)]

t1 <- t[1:tbreak]
t2 <- t[tbreak:length(t)]

o1<-GetRho(x1,t1)
o2<-GetRho(x2,t2)

mu1 <- mean(x1)
sigma1 <- sd(x1)
rho1 <- o1[1]

mu2 <- mean(x2)
sigma2 <- sd(x2)
rho2 <- o2[1]

LL1 <- -o1[2]
LL2 <- -o2[2]

m <- matrix(c(mu1,mu2,sigma1,sigma2,rho1,rho2,LL1,LL2),2,4)
colnames(m) <- c("mu","sigma","rho","LL")

return(m)
}

SweepBreaks <- function(x,t,range=0.6)
{
n<-length(t)
start <- (1-range)/2
breaks<-round((start*n):((1-start)*n))
Ls<-breaks*0

l<-length(breaks)

BreakParams <- matrix(NA,l,8)
#BreakParams <- data.frame(mu1=NA,s1=NA,rho1=NA,LL1=NA,mu2=NA,s2=NA,rho2=NA,LL2=NA)

for(i in 1:l)
{
myDoubleL <- GetDoubleL(x,t,breaks[i])
BreakParams[i,] <- c(myDoubleL[1,],myDoubleL[2,])
}

# remember: LL1 and LL2 are columns 4 and 8
BreakMatrix<- cbind(breaks,t[breaks], BreakParams,
BreakParams[,4]+BreakParams[,8])

colnames(BreakMatrix) <- c("breaks","tbreaks","mu1","sigma1","rho1","LL1","mu2","sigma2","rho2","LL2","LL")

return(BreakMatrix[2:nrow(BreakMatrix),])
}

GetModels <- function(x,t,tbreak,K=2)
{
for(i in 0:7)
{
f<-get(paste("M",i,sep=""))
myr<-data.frame(Model=i,f(x,t,tbreak,K))
ifelse(i==0,
r<-myr,
r<-rbind(r,myr))
}
return(r)
}

M0 <- function(x,t,tbreak,K=2)
# null model: all mus, s's, rhos the same
{
rhoLL <- GetRho(x,t)
LL <- rhoLL[2]

bic <- K*LL + 3*log(length(x))

rho1<-rhoLL[1]
rho2<-rho1

mu1<-mean(x)
mu2<-mu1

s1<-sd(x)
s2<-s1

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M1 <- function(x,t,tbreak,K=2)
# mus different, all else the same
{

x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x1)
mu2<-mean(x2)

xprime <- c(x1-mu1,x2-mu2)
s1<-sd(xprime)
s2<-s1

rho1<-as.numeric(GetRho(xprime,t)[1])
rho2<-rho1

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)
LL<-LL1+LL2
bic <- K*LL + 5*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M2 <- function(x,t,tbreak,K=2)
# sds different, all else same
{

x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]

t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x)
mu2<-mu1

s1<-sd(x1)
s2<-sd(x2)

xprime <- c( (x1-mu1)/s1 , (x2-mu2)/s2 )
rho1 <- as.numeric(GetRho(xprime,t)[1])
rho2 <- rho1

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)
LL<-LL1+LL2
bic <- K*LL + 5*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M3 <- function(x,t,tbreak,K=2)
# rhos different, all else same
{
x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1 <- mean(x)
mu2 <- mu1

s1 <- sd(x)
s2 <- s1

rho1<-as.numeric(GetRho(x1,t1)[1])
rho2<-as.numeric(GetRho(x2,t2)[1])

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)
LL<-LL1+LL2
bic <- K*LL + 5*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M4 <- function(x,t,tbreak,K=2)
# mu and sigma different, rho same
{
x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x1)
mu2<-mean(x2)
s1<-sd(x1)
s2<-sd(x2)

xprime <- c( (x1-mu1)/s1 , (x2-mu2)/s2 )
rho1 <- as.numeric(GetRho(xprime,t)[1])
rho2 <- rho1

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)
LL<-LL1+LL2
bic <- K*LL + 6*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M5 <- function(x,t,tbreak,K=2)
# mu and rho different, sigma same
{
x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x1)
mu2<-mean(x2)

xprime <- c(x1-mu1, x2-mu2)
s1<-sd(xprime)
s2<-s1

rho1<-as.numeric(GetRho(x1,t1)[1])
rho2<-as.numeric(GetRho(x2,t2)[1])

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)
LL<-LL1+LL2
bic <- K*LL+ 6*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M6 <- function(x,t,tbreak,K=2)
# sigma and rho different, mu same
{
x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x)
mu2<-mean(x)
s1<-sd(x1)
s2<-sd(x2)

x1prime <- (x1-mu1)/s1
x2prime <- (x2-mu2)/s2

rho1 <- as.numeric(GetRho(x1prime,t1)[1])
rho2 <- as.numeric(GetRho(x2prime,t2)[1])

LL1<-GetLL(x1,t1,mu1,s1,rho1)
LL2<-GetLL(x2,t2,mu2,s2,rho2)

LL<-LL1+LL2

bic <- K*LL+ 6*log(length(x))

return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

M7 <- function(x,t,tbreak,K=2)
# most "alternative" model: all mus, s's, rhos different
{
rhoLL1 <- GetRho(x[1:tbreak],t[1:tbreak])
rhoLL2 <- GetRho(x[(tbreak+1):length(x)],t[(tbreak+1):length(x)])

LL1 <- rhoLL1[2]
LL2 <- rhoLL2[2]

x1<-x[1:tbreak]
x2<-x[(tbreak+1):length(x)]
t1<-t[1:tbreak]
t2<-t[(tbreak+1):length(x)]

mu1<-mean(x1)
mu2<-mean(x2)
s1<-sd(x1)
s2<-sd(x2)
rho1 <- rhoLL1[1]
rho2 <- rhoLL2[1]

LL <- LL1+LL2
bic <- K*LL + 7*log(length(x))
return(data.frame(LL,bic,mu1,s1,rho1,mu2,s2,rho2))
}

WindowSweep <- function (x, t, windowsize = 50, windowstep = 1, sine = 0, K = 2, plotme = TRUE)
{
low <- seq(1, (length(t) - windowsize), windowstep)
hi <- low + windowsize
for (i in 1:length(low)) {
myx <- x[low[i]:hi[i]]
myt <- t[low[i]:hi[i]]
bp <- SweepBreaks(myx, myt)
myestimate <- bp[bp[, 11] == max(bp[, 11]), ]
breakpoint <- myestimate[1]
tbreak <- myestimate[2]
ifelse(sine,
allmodels <- GetModelsSin(myx, myt, breakpoint, K),
allmodels <- GetModels(myx, myt, breakpoint, K))
mymodel <- allmodels[allmodels\$bic == min(allmodels\$bic),]
mymodel <- data.frame(mymodel, Break = tbreak)
ifelse(i == 1,
estimates <- mymodel,
estimates <- rbind(estimates, mymodel))
if (plotme) {
plot.ts(t, x, type = "l", col = "grey")
lines(t, x, type = "l")
lines(myt, myx, col = "green")
abline(v = tbreak)
print(estimates[i, ])
}
}
return(data.frame(estimates))
}

PartitionParameters <- function(ws,t,windowsize=50,windowstep=1)
{
n.col<-length(t)
n.row<-dim(ws)[1]

mu.M <- matrix(NA,n.row,n.col)
s.M <- matrix(NA,n.row,n.col)
rho.M <- matrix(NA,n.row,n.col)

for(i in 1:n.row)
{
myws<-ws[i,]
dts <- abs(t-myws\$Break)
tbreak <- match(min(dts),dts)

max <- min(n.col,i+windowsize)

mu.M[i,i:tbreak] <- myws\$mu1
mu.M[i,(tbreak+1):max] <- myws\$mu2
s.M[i,i:tbreak] <- myws\$s1
s.M[i,(tbreak+1):max] <- myws\$s2
rho.M[i,i:tbreak] <- myws\$rho1
rho.M[i,(tbreak+1):max] <- myws\$rho2
}

return(data.frame(mu.hat,s.hat,rho.hat))
}

PlotBCPA <- function (t, x, ws, pp, threshold=10)
{
plot(t, x, type = "n", main = "", ylim = c(min(x), max(x) *
1.3))
x.breaks <- ws\$Break
x.model <- ws\$Model
x.breaks <- x.breaks[x.model > 0]
x.model <- x.model[x.model > 0]
mids <- hist(x.breaks, breaks = t, plot = F)\$mid
freq <- hist(x.breaks, breaks = t, plot = F)\$count
goodbreaks <- mids[freq > threshold]
freq <- freq[freq > threshold]
break.cols <- heat.colors(max(freq))
break.cols <- break.cols[length(break.cols):1]
abline(v = goodbreaks, lwd = freq, col = "orange")
lines(t, x, col = "darkgrey")
lines(t, pp\$mu.hat, lwd = 2)
lines(t, pp\$mu.hat + pp\$s.hat, col = 2, lwd = 1.5)
lines(t, pp\$mu.hat - pp\$s.hat, col = 2, lwd = 1.5)
rho.hat <- pp\$rho.hat
rho.int <- round((rho.hat - min(rho.hat, na.rm = 1)) * 1000) +
1
rho.col <- topo.colors(max(rho.int, na.rm=TRUE))[rho.int]
rho.cex <- (rho.hat - min(rho.hat, na.rm=TRUE))/(diff(range(rho.hat))) *
1.5 + 0.5
points(t, x, col = rho.col, pch = 19)
rho.index <- quantile(0:max(rho.int))
rho.index[1] <- 1
legend.cols <- topo.colors(max(rho.int, na.rm=TRUE))[round(rho.index)]
legend.rhos <- round(quantile((min(rho.hat) * 1000):(max(rho.hat) *
1000))/1000, 2)
legend("bottomright", bg = "white", legend = c(expression(hat(rho)),
legend.rhos), pch = c(0, rep(19, 5)), ncol = 2, col = c(0,
legend.cols), xjust = 0.5, yjust = 0.5, cex = 1.2)
legend("topright", bg = "white", legend = c(expression(hat(mu)),
expression(hat(mu) %+-% hat(sigma))), lty = 1, lwd = 2:1,
col = 1:2, xjust = 0.5, yjust = 0.5, cex = 1.2)
}

PathPlot <- function(t, z, ws, pp, threshold = 10)
{
n<-names(pp)
for(i in 1:length(n))
assign(n[i],pp[,i])

# Prep rho legend
rho.int <- round((rho.hat-min(rho.hat,na.rm=1))*1000)+1
rho.col <- topo.colors(max(rho.int))[rho.int]
rho.index<-quantile(0:max(rho.int))
rho.index[1]<-1

legend.cols<-topo.colors(max(rho.int))[round(rho.index)]
legend.rhos<-round(quantile((min(rho.hat)*1000):(max(rho.hat)*1000))/1000,1)

# Prep mu legend
sizes <- round(quantile(seq(0,max(mu.hat,na.rm=1),.01)),2)[2:5]

# get "notable" changepoints and corresponding locations
x.breaks <- ws\$Break
x.model <- ws\$Model

mids <- hist(x.breaks, breaks = t, plot = F)\$mid
freq <- hist(x.breaks, breaks = t, plot = F)\$count

t.breaks <- mids[freq > threshold]
z.mid <- (z[-(1:2)]+z[-( (length(z)-1):(length(z)))])/2
t.mid <- (t[-1]+t[-length(t)])/2
z.breaks <- z.mid[match(t.breaks, t.mid)+1]

# plot everything
plot(z, asp=1, type="l", xlab="X", ylab="Y")
segments(Re(z[-1]), Im(z[-1]), Re(z[-length(z)]), Im(Z[-length(Z)]),
lwd=mu.hat/max(mu.hat)*12, col="darkgrey")
segments(Re(z[-1]), Im(z[-1]), Re(z[-length(z)]), Im(Z[-length(Z)]),
lwd=mu.hat/max(mu.hat)*10, col=rho.col)
points(z.breaks, col=rgb(1,.2,.2,.8), cex=sqrt(freq/max(freq))*5, pch=4, lwd=2)
points(z,pch=19,cex=0.5)
lines(z, lwd=0.75)