BCPA/All BCPA Functions
From QERM Wiki
(Difference between revisions)
| Line 413: | Line 413: | ||
PathPlot <- function(t, z, ws, pp, threshold = 10) | PathPlot <- function(t, z, ws, pp, threshold = 10) | ||
{ | { | ||
| − | + | n<-names(pp) | |
| − | + | for(i in 1:length(n)) | |
| − | assign(n[i], | + | assign(n[i],pp[,i]) |
# Prep rho legend | # Prep rho legend | ||
| Line 448: | Line 448: | ||
segments(Re(z[-1]), Im(z[-1]), Re(z[-length(z)]), Im(Z[-length(Z)]), | 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) | lwd=mu.hat/max(mu.hat)*10, col=rho.col) | ||
| − | points(z.breaks, col=rgb(1,.2,.2,.8), cex=sqrt( | + | 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) | points(z,pch=19,cex=0.5) | ||
lines(z, lwd=0.75) | lines(z, lwd=0.75) | ||
| Line 458: | Line 458: | ||
lty=1, lwd=sizes/max(mu.hat)*10,title=expression(hat(mu)),cex=1.2,bty="n") | lty=1, lwd=sizes/max(mu.hat)*10,title=expression(hat(mu)),cex=1.2,bty="n") | ||
} | } | ||
| − | |||
</pre> | </pre> | ||
[[Category: BCPA]] | [[Category: BCPA]] | ||
Latest revision as of 12:15, 17 August 2012
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
}
adjust <- colSums(!is.na(mu.M))
mu.hat<-colSums(mu.M,na.rm=1)/adjust
s.hat<-colSums(s.M,na.rm=1)/adjust
rho.hat<-colSums(rho.M,na.rm=1)/adjust
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)
# add a legend
legend("topright",fill=legend.cols,legend=legend.rhos,
title=expression(hat(rho)),bty="n",cex=1.2)
legend("bottomleft",col="darkgrey",legend=round(sizes, -1),
lty=1, lwd=sizes/max(mu.hat)*10,title=expression(hat(mu)),cex=1.2,bty="n")
}
