BCPA/All BCPA Functions

From QERM Wiki
Jump to: navigation, search
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")
}
Personal tools