Nothing
R/plot.mrbeast.R
In Rbeast: Bayesian Change-Point Detection and Time Series Decomposition
Defines functions
plot.mrbeast
plot.mrbeast<-function(
x,
index = 1,
vars = c('st','s','scp','sorder','t','tcp','torder','o','ocp','error'),
col = NULL,
main ="BEAST decomposition and changepoint detection",
xlab ='Time',
ylab = NULL,
cex.main =1,
cex.lab =1,
relative.heights= NULL,
interactive = FALSE,
ncpStat = c('median','mode','mean','pct90','max'),
... )
{
ncpStat=match.arg(ncpStat)
#vars = c('st','s','scp','sorder','t','tcp','torder','error')
#vars = c('st','s','scp','sorder','samp','t','tcp','torder', 'tslp','o','ocp','error'),
#col = NULL
#main ="BEAST decomposition and changepoint detection"
#xlab ='Time'
#ylab = NULL
#cex.main = 1.2
#cex.lab = 1
#relative.heights=NULL
#interactive=FALSE
vars = tolower(vars)
vars_log = vars=='st'|vars=='s'|vars=='t'|vars=='scp'|vars=='tcp'|vars=='sorder'|vars=='torder'|vars=='error'|vars=='o'|vars=='ocp'|vars=='samp'|vars=='tslp'|vars=='slpsgn';
vars = vars[vars_log]
if (length(ylab)== length(vars_log)){
ylab = ylab[vars_log]
} else{
ylab = vector(mode='character',length(vars))
ylab[vars=='st'] ='Y'
ylab[vars=='s'] ='season'
ylab[vars=='t'] ='trend'
ylab[vars=='o'] ='outlier'
ylab[vars=='scp'] ="Pr(scp)"
ylab[vars=='tcp'] ="Pr(tcp)"
ylab[vars=='ocp'] ="Pr(ocp)"
ylab[vars=='sorder']= expression("order"[S])
ylab[vars=='torder']= expression("order"[T])
ylab[vars=='samp'] = 'amplitude'
ylab[vars=='tslp'] = 'slope'
ylab[vars=='slpsgn'] = "slpSign"
ylab[vars=='error'] = "error"
}
if (length(col)== length(vars_log)){
col = col[vars_log]
} else{
col = vector(mode='character',length(vars))
col[vars=='st'] ='#111111'
col[vars=='s'] ='red'
col[vars=='scp'] ='red'
col[vars=='sorder']='red'
col[vars=='samp'] ='red'
col[vars=='t'] ='green'
col[vars=='tcp'] ='green'
col[vars=='torder']='green'
col[vars=='tslp'] ='green'
col[vars=='slpsgn'] ='black' #the border only
col[vars=='o'] ='blue'
col[vars=='ocp'] ='blue'
col[vars=='error'] ='darkgray'
}
if ( length(x$marg_lik)> 1 ) {
x=tsextract(x,index) # more than 1 time series is present
}
q = length(x$trend$Y)
heights=relative.heights
if (length(heights)== length(vars_log)){
heights = heights[vars_log]
} else{
heights = (1:length(vars))
heights =heights-heights+1
heights[vars=='st']=.8/q
heights[vars=='s']= .8/q
heights[vars=='t']=.8/q
heights[vars=='o']=.8/q
heights[vars=='scp']=.2
heights[vars=='tcp']=.2
heights[vars=='ocp']=.2
heights[vars=='sorder']=.2
heights[vars=='torder']=.2
heights[vars=='samp'] =.4/q
heights[vars=='tslp'] =.4/q
heights[vars=='slpsgn'] =.4/q
heights[vars=='error'] =0.4/q
}
hasSeason =!is.null(x$season)
hasOutlier =!is.null(x$outlier)
hasHarmonic =!is.null(x$season$order)
hasTOrder =!is.null(x$trend$order)
hasData =!is.null(x$data)
hasAmp =!is.null(x$season$amp)
hasSlp =!is.null(x$trend$slp)
idx =1:length(vars) > 0
if(!hasAmp) { idx = idx & !(vars=='samp')}
if(!hasSlp) { idx = idx & !(vars=='tslp'|vars=='slpsgn')}
if(!hasSeason) { idx = idx & !(vars=='st'|vars=='s'|vars=='sorder'|vars=='scp') }
if(!hasHarmonic) { idx = idx & !(vars=='sorder') }
if(!hasTOrder) { idx = idx & !(vars=='torder') }
if(!hasOutlier) { idx = idx &!(vars=='o'|vars=='ocp') }
if(!hasData) { idx = idx &!(vars=='error') }
col = col[idx]
vars = vars[idx]
ylab = ylab[idx]
heights = heights[idx]
nPlots=length(vars)
newHeights=NULL;
for (i in 1:nPlots){
v=vars[i]
if (v=='st' || v=='s' || v== 't' || v == 'o' || v == 't' || v == 'tslp' || v == 'samp' || v == 'error') {
newHeights=c(newHeights, rep( heights[i], q) );
} else{
newHeights=c(newHeights, heights[i]);
}
}
heights=newHeights;
newYlab=NULL;
for (i in 1:nPlots){
v=vars[i]
if (v=='st' || v=='s' || v== 't' || v == 'o' || v == 't' || v == 'tslp' || v == 'samp' || v == 'error') {
newYlab=c(newYlab, rep( ylab[i], q) );
} else{
newYlab=c(newYlab, ylab[i]);
}
}
ylab=newYlab;
if(nPlots==0){
stop("No valid variable names speciffied int the 'vars' argument. Possible names include 'st','t','s','sorder','torder','scp','tcp','samp','tslp','o', 'ocp', and 'error'. ");
}
#######################################################
# Functions and variables to load the outputs
########################################################
t = x$time;
data = x$data
t2t = c(t, rev(t))
N = length(t)
CI = 0
Y = 0
SD = 0
Prob2Prob = 0;
Prob = 0;
Order = 0
Slp =0
SlpSD =0
SlpSignPos =0
SlpSignZero=0
cp =0
cpCI =0
ncp =0
ncpPr=0
cpPr =0
cpChange=0
HasChangePoint=FALSE;
":" = function(i,j){ if(i>j){return(NULL)} ; seq(i,j) }
Yts = 0
YtsSD = 0
Yerr = 0
get.Yts = function() {
sig2 = diag(x$sig2)
Yts <<- list();
YtsSD <<- list();
Yerr <<- list();
for (i in 1:q) {
Yts[[i]] <<- x$trend$Y[[i]]
SD2 = x$trend$SD[[i]]^2 +sig2[i]
if (!is.null(x$season)) {
Yts[[i]] <<- Yts[[i]]+x$season$Y[[i]]
SD2 = SD2+x$season$SD[[i]]^2
}
if (!is.null(x$outlier)) {
Yts[[i]] <<- Yts[[i]]+x$outlier$Y[[i]]
SD2 = SD2+x$outlier$SD[[i]]^2
}
SD = sqrt(SD2)
YtsSD[[i]] <<- c(Yts[[i]]-SD, rev(Yts[[i]]+SD) )
Yerr[[i]] <<- data[[i]]-Yts[[i]]
}
}
get.T = function(){
Y <<- list()
SD <<- list()
CI <<- list()
Slp <<- list()
SlpSD <<- list()
SlpSignPos <<- list()
SlpSignZero <<- list()
for (i in 1:q){
Y[[i]] <<-x$trend$Y[[i]];
SD[[i]] <<-c(Y[[i]]-x$trend$SD[[i]], rev(Y[[i]]+x$trend$SD[[i]]));
if ( !is.null(x$trend$CI) ) CI[[i]] <<-c(x$trend$CI[[i]][,1], rev(x$trend$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
if (hasSlp){
Slp[[i]] <<- x$trend$slp[[i]]
SlpSD[[i]] <<- c(Slp[[i]]-x$trend$slpSD[[i]], rev(Slp[[i]]+x$trend$slpSD[[i]]));
SlpSignPos[[i]] <<- x$trend$slpSgnPosPr[[i]]
SlpSignZero[[i]] <<- x$trend$slpSgnZeroPr[[i]]
}
}
Order <<- x$trend$order;
HasChangePoint <<- !is.null(x$trend$cp);
}
get.tcp = function(){
cmpnt = x$trend
cp <<- cmpnt$cp;
cpCI <<- cmpnt$cpCI;
ncp <<- switch(ncpStat, mode=cmpnt$ncp_mode, median=cmpnt$ncp_median,mean=cmpnt$ncp,pct90=cmpnt$ncp_pct90,max=sum(!is.nan(cp)))
ncp <<- base::round(ncp)
ncpPr <<- cmpnt$ncpPr
cpPr <<- cmpnt$cpPr
cpChange<<- cmpnt$cpAbruptChange
HasChangePoint <<- !is.null(x$trend$cp);
Prob <<- cmpnt$cpOccPr;
Prob2Prob <<- c(Prob,Prob-Prob)
}
###########################################################
Amp = 0
AmpSD = 0
get.S = function(){
Y <<- list()
SD <<- list()
CI <<- list()
Amp <<- list()
AmpSD <<- list()
for (i in 1:q){
Y[[i]] <<-x$season$Y[[i]];
SD[[i]] <<-c(Y[[i]] -x$season$SD[[i]], rev(Y[[i]] +x$season$SD[[i]]));
if ( !is.null(x$season$CI) ) CI[[i]] <<-c(x$season$CI[[i]][,1], rev(x$season$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
if (hasAmp){
Amp[[i]] <<-x$season$amp[[i]]
AmpSD[[i]] <<-c(Amp[[i]]-x$season$ampSD[[i]], rev(Amp[[i]]+x$season$ampSD[[i]]));
}
}
Order <<-x$season$order;
HasChangePoint <<- !is.null(x$season$cp);
}
get.scp = function(){
cmpnt = x$season
cp <<-cmpnt$cp;
cpCI <<-cmpnt$cpCI;
ncp <<-switch(ncpStat, mode=cmpnt$ncp_mode, median=cmpnt$ncp_median,mean=cmpnt$ncp,pct90=cmpnt$ncp_pct90,max=sum(!is.nan(cp)))
ncp <<-base::round(ncp)
ncpPr <<-cmpnt$ncpPr
cpPr <<-cmpnt$cpPr
cpChange<<-cmpnt$cpAbruptChange
Prob <<-x$season$cpOccPr;
Prob2Prob <<-c(Prob,Prob-Prob)
HasChangePoint <<- !is.null(x$season$cp);
}
get.O = function(){
Y <<-list()
SD <<-list()
for (i in 1:q){
Y[[i]] <<-x$outlier$Y[[i]] ;
SD[[i]] <<-c(Y[[i]] -x$outlier$SD[[i]] , rev(Y[[i]] +x$outlier$SD[[i]] ));
if ( !is.null(x$outlier$CI) ) CI[[i]] <<-c(x$outlier$CI[[i]][,1], rev(x$outlier$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
}
}
get.ocp = function(){
#cp <<-x$season$cp;
#cpCI <<-x$season$cpCI;
#ncp <<-sum(!is.nan(cp))
#ncpPr<<-x$season$ncpPr
#cpPr <<-x$season$cpPr
Prob <<-x$outlier$cpOccPr;
Prob2Prob <<-c(Prob,Prob-Prob)
}
YSIDE = 2;
padj = -1.2
tcl = 0.2
plotId = 1;
plot.st=function(col=c(0.2,0.2,0.2), ylabel){
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
if (hasData){
Yall=c(YtsSD[[i]],data[[i]])
plot( c(t[1],t[N]), c(min(Yall,na.rm = TRUE),max(Yall,na.rm = TRUE)), type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
} else {
plot( t2t, YtsSD[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
}
### polygon(t2t, YtsSD[[i]], col = rgb(.5,0.5,0.5,.5), border = NA);
if (hasData){
points( t, data[[i]], type = 'p', col='#777777');
}
points( t, Yts[[i]], type = 'l', lwd=2,col=rgb(col[1],col[2],col[3]));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.y =function(col=c(0,1,0), ylabel){
alpha=0.1
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
if (hasData && !hasSeason){
Yall=c(CI[[i]],data[[i]])
plot( c(t[1],t[N]), c(min(Yall,na.rm = TRUE),max(Yall,na.rm = TRUE)), type = 'n',ann=FALSE, xaxt='n', yaxt='n');
} else {
plot( t2t, CI[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
}
if (hasData && !hasSeason){
points(t, data[[i]], type = 'p', col='#777777');
}
### polygon(t2t, CI[[i]], col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Y[[i]], type = 'l',col=rgb(col[1],col[2],col[3] ) );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.prob <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
alpha=0.2
plot( c(t2t[1],t2t), c(0.22,Prob2Prob),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
polygon(t2t, Prob2Prob, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Prob, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.order <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
maxOrder=max(1.01,max(Order));
plot( c(t[1],t[length(t)]), c(0,maxOrder),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
points( t, Order, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.amp =function(col=c(0,1,0), ylabel){
alpha=0.5
plot( t, Amp, type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, AmpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Amp, type = 'l',col=rgb(col[1],col[2],col[3] ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.slp =function(col=c(0,1,0), ylabel){
alpha=0.5
plot( t2t, SlpSD, type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
polygon(t2t, SlpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Slp, type = 'l',col=rgb(0,0,0,0.8 ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.slpsgn = function(col=c(0,1,0), ylabel){
alpha=0.5
SlpSignNeg = 1-SlpSignPos-SlpSignZero
plot( c(t[1],t[length(t)]), c(0,1), type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, SlpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
y2y = c(t-t, rev(SlpSignNeg)); polygon(t2t, y2y, col=rgb(0,0,1,alpha),border=NA)
y2y = c(SlpSignNeg, rev(1-SlpSignPos)); polygon(t2t, y2y, col=rgb(0,1,0,alpha),border=NA)
y2y = c(1-SlpSignPos, t-t+1); polygon(t2t, y2y, col=rgb(1,0,0,alpha),border=NA)
points( t, t-t+0.5 , type = 'l', lty=2, col=rgb(col[1],col[2],col[3] ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.o =function(col=c(0,1,0), ylabel){
alpha=0.5
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
plot( t2t, CI[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, CI, col = rgb(col[1],col[2],col[3],alpha), border = NA);
#points( t, Y, type = 'l',col='#333333');
segments(t,t-t,t,Y[[i]], col= rgb(col[1],col[2],col[3]))
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.oprob <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
alpha=0.2
plot( c(t2t[1],t2t), c(0.22,Prob2Prob),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
# polygon(t2t, Prob2Prob, col = rgb(col[1],col[2],col[3],alpha), border = NA);
# points( t, Prob, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
segments(t,t-t,t,Prob, col= rgb(col[1],col[2],col[3]),lwd=1.5)
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.error=function(col=c(0.2,0.2,0.2), ylabel){
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
plot( t, Yerr[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
lines(t, t-t, col= rgb(col[1],col[2],col[3]))
segments(t,t-t,t,Yerr[[i]], col= rgb(col[1],col[2],col[3]))
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.new()
opar= par(c('oma','mar','fig'))
par(oma=c(0,0,0,0))
par(mar=c(0,2,0,2))
par(fig=c(0,1,0,1))
lm= strheight('a',units='figure')*dev.size()[2]/dev.size()[1]*4.5
rm= strheight('a',units='figure')*dev.size()[2]/dev.size()[1]*2.3
bm= strheight('a',units='figure')*4.0
tm= strheight('a',units='figure')*4.0
#######################################################
# Create a subplot given the relative heights of the vertical plots
########################################################
subplot <- function(i, heights){
Hmiddle =1-tm-bm
h = (heights/sum(heights))*Hmiddle
hcum = cumsum(h)
plt = c(lm, 1-rm, 1-tm-hcum[i], 1-tm-hcum[i]+h[i])
par(plt=plt,new=TRUE)
}
for (i in 1:length(vars)){
ytitle = ylab[i]
var = tolower(vars[i])
clr = as.numeric(col2rgb( col[i] ))/255
if (var=='st') { get.Yts(); plot.st(clr, ytitle) }
if (var=='s' ) { get.S(); get.scp(); plot.y(clr, ytitle) }
if (var=='t' ) { get.T(); get.tcp(); plot.y(clr, ytitle) }
if (var=='scp') { get.S(); get.scp(); plot.prob(clr, ytitle) }
if (var=='tcp') { get.T(); get.tcp(); plot.prob(clr, ytitle) }
if (var=='sorder'){ get.S(); get.scp(); plot.order(clr, ytitle) }
if (var=='torder'){ get.T(); get.tcp(); plot.order(clr, ytitle) }
if (var=='samp') { get.S(); plot.amp(clr, ytitle) }
if (var=='tslp') { get.T(); plot.slp(clr, ytitle) }
if (var=='slpsgn'){ get.T(); plot.slpsgn(clr, ytitle) }
if (var=='o') { get.O(); plot.o(clr, ytitle) }
if (var=='ocp') { get.ocp(); plot.oprob(clr, ytitle) }
if (var=='error') { get.Yts(); plot.error(clr, ytitle) }
#title(ylab = ytitle,cex=cex.lab, mgp = c(1.25, 1.25, 0));
if(i==1){
#mtext(main,line=.5,cex=cex.main,font=2)
}
if(i==nPlots){
title(xlab = xlab ,mgp = c(1.2, 1.2, 0));
if (sum(class(t)=='Date') >0 ){
axis.Date(1, t, labels = TRUE, padj =padj, tcl =tcl,cex.axis=1);
} else{
axis( 1, labels = TRUE, padj =padj, tcl =tcl,cex.axis=1);
}
}
}
par(opar)
#text(mean(time) / 2.0, 'NO SEASONAL COMPONET.IGNORE THIS POLOT');
}
#plot.mrbeast(o,vars=c('st','s','scp','t','tcp','o','ocp'))
Try the Rbeast package in your browser
Any scripts or data that you put into this service are public.
Rbeast documentation built on Sept. 12, 2024, 7:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.mrbeast
plot.mrbeast<-function(
x,
index = 1,
vars = c('st','s','scp','sorder','t','tcp','torder','o','ocp','error'),
col = NULL,
main ="BEAST decomposition and changepoint detection",
xlab ='Time',
ylab = NULL,
cex.main =1,
cex.lab =1,
relative.heights= NULL,
interactive = FALSE,
ncpStat = c('median','mode','mean','pct90','max'),
... )
{
ncpStat=match.arg(ncpStat)
#vars = c('st','s','scp','sorder','t','tcp','torder','error')
#vars = c('st','s','scp','sorder','samp','t','tcp','torder', 'tslp','o','ocp','error'),
#col = NULL
#main ="BEAST decomposition and changepoint detection"
#xlab ='Time'
#ylab = NULL
#cex.main = 1.2
#cex.lab = 1
#relative.heights=NULL
#interactive=FALSE
vars = tolower(vars)
vars_log = vars=='st'|vars=='s'|vars=='t'|vars=='scp'|vars=='tcp'|vars=='sorder'|vars=='torder'|vars=='error'|vars=='o'|vars=='ocp'|vars=='samp'|vars=='tslp'|vars=='slpsgn';
vars = vars[vars_log]
if (length(ylab)== length(vars_log)){
ylab = ylab[vars_log]
} else{
ylab = vector(mode='character',length(vars))
ylab[vars=='st'] ='Y'
ylab[vars=='s'] ='season'
ylab[vars=='t'] ='trend'
ylab[vars=='o'] ='outlier'
ylab[vars=='scp'] ="Pr(scp)"
ylab[vars=='tcp'] ="Pr(tcp)"
ylab[vars=='ocp'] ="Pr(ocp)"
ylab[vars=='sorder']= expression("order"[S])
ylab[vars=='torder']= expression("order"[T])
ylab[vars=='samp'] = 'amplitude'
ylab[vars=='tslp'] = 'slope'
ylab[vars=='slpsgn'] = "slpSign"
ylab[vars=='error'] = "error"
}
if (length(col)== length(vars_log)){
col = col[vars_log]
} else{
col = vector(mode='character',length(vars))
col[vars=='st'] ='#111111'
col[vars=='s'] ='red'
col[vars=='scp'] ='red'
col[vars=='sorder']='red'
col[vars=='samp'] ='red'
col[vars=='t'] ='green'
col[vars=='tcp'] ='green'
col[vars=='torder']='green'
col[vars=='tslp'] ='green'
col[vars=='slpsgn'] ='black' #the border only
col[vars=='o'] ='blue'
col[vars=='ocp'] ='blue'
col[vars=='error'] ='darkgray'
}
if ( length(x$marg_lik)> 1 ) {
x=tsextract(x,index) # more than 1 time series is present
}
q = length(x$trend$Y)
heights=relative.heights
if (length(heights)== length(vars_log)){
heights = heights[vars_log]
} else{
heights = (1:length(vars))
heights =heights-heights+1
heights[vars=='st']=.8/q
heights[vars=='s']= .8/q
heights[vars=='t']=.8/q
heights[vars=='o']=.8/q
heights[vars=='scp']=.2
heights[vars=='tcp']=.2
heights[vars=='ocp']=.2
heights[vars=='sorder']=.2
heights[vars=='torder']=.2
heights[vars=='samp'] =.4/q
heights[vars=='tslp'] =.4/q
heights[vars=='slpsgn'] =.4/q
heights[vars=='error'] =0.4/q
}
hasSeason =!is.null(x$season)
hasOutlier =!is.null(x$outlier)
hasHarmonic =!is.null(x$season$order)
hasTOrder =!is.null(x$trend$order)
hasData =!is.null(x$data)
hasAmp =!is.null(x$season$amp)
hasSlp =!is.null(x$trend$slp)
idx =1:length(vars) > 0
if(!hasAmp) { idx = idx & !(vars=='samp')}
if(!hasSlp) { idx = idx & !(vars=='tslp'|vars=='slpsgn')}
if(!hasSeason) { idx = idx & !(vars=='st'|vars=='s'|vars=='sorder'|vars=='scp') }
if(!hasHarmonic) { idx = idx & !(vars=='sorder') }
if(!hasTOrder) { idx = idx & !(vars=='torder') }
if(!hasOutlier) { idx = idx &!(vars=='o'|vars=='ocp') }
if(!hasData) { idx = idx &!(vars=='error') }
col = col[idx]
vars = vars[idx]
ylab = ylab[idx]
heights = heights[idx]
nPlots=length(vars)
newHeights=NULL;
for (i in 1:nPlots){
v=vars[i]
if (v=='st' || v=='s' || v== 't' || v == 'o' || v == 't' || v == 'tslp' || v == 'samp' || v == 'error') {
newHeights=c(newHeights, rep( heights[i], q) );
} else{
newHeights=c(newHeights, heights[i]);
}
}
heights=newHeights;
newYlab=NULL;
for (i in 1:nPlots){
v=vars[i]
if (v=='st' || v=='s' || v== 't' || v == 'o' || v == 't' || v == 'tslp' || v == 'samp' || v == 'error') {
newYlab=c(newYlab, rep( ylab[i], q) );
} else{
newYlab=c(newYlab, ylab[i]);
}
}
ylab=newYlab;
if(nPlots==0){
stop("No valid variable names speciffied int the 'vars' argument. Possible names include 'st','t','s','sorder','torder','scp','tcp','samp','tslp','o', 'ocp', and 'error'. ");
}
#######################################################
# Functions and variables to load the outputs
########################################################
t = x$time;
data = x$data
t2t = c(t, rev(t))
N = length(t)
CI = 0
Y = 0
SD = 0
Prob2Prob = 0;
Prob = 0;
Order = 0
Slp =0
SlpSD =0
SlpSignPos =0
SlpSignZero=0
cp =0
cpCI =0
ncp =0
ncpPr=0
cpPr =0
cpChange=0
HasChangePoint=FALSE;
":" = function(i,j){ if(i>j){return(NULL)} ; seq(i,j) }
Yts = 0
YtsSD = 0
Yerr = 0
get.Yts = function() {
sig2 = diag(x$sig2)
Yts <<- list();
YtsSD <<- list();
Yerr <<- list();
for (i in 1:q) {
Yts[[i]] <<- x$trend$Y[[i]]
SD2 = x$trend$SD[[i]]^2 +sig2[i]
if (!is.null(x$season)) {
Yts[[i]] <<- Yts[[i]]+x$season$Y[[i]]
SD2 = SD2+x$season$SD[[i]]^2
}
if (!is.null(x$outlier)) {
Yts[[i]] <<- Yts[[i]]+x$outlier$Y[[i]]
SD2 = SD2+x$outlier$SD[[i]]^2
}
SD = sqrt(SD2)
YtsSD[[i]] <<- c(Yts[[i]]-SD, rev(Yts[[i]]+SD) )
Yerr[[i]] <<- data[[i]]-Yts[[i]]
}
}
get.T = function(){
Y <<- list()
SD <<- list()
CI <<- list()
Slp <<- list()
SlpSD <<- list()
SlpSignPos <<- list()
SlpSignZero <<- list()
for (i in 1:q){
Y[[i]] <<-x$trend$Y[[i]];
SD[[i]] <<-c(Y[[i]]-x$trend$SD[[i]], rev(Y[[i]]+x$trend$SD[[i]]));
if ( !is.null(x$trend$CI) ) CI[[i]] <<-c(x$trend$CI[[i]][,1], rev(x$trend$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
if (hasSlp){
Slp[[i]] <<- x$trend$slp[[i]]
SlpSD[[i]] <<- c(Slp[[i]]-x$trend$slpSD[[i]], rev(Slp[[i]]+x$trend$slpSD[[i]]));
SlpSignPos[[i]] <<- x$trend$slpSgnPosPr[[i]]
SlpSignZero[[i]] <<- x$trend$slpSgnZeroPr[[i]]
}
}
Order <<- x$trend$order;
HasChangePoint <<- !is.null(x$trend$cp);
}
get.tcp = function(){
cmpnt = x$trend
cp <<- cmpnt$cp;
cpCI <<- cmpnt$cpCI;
ncp <<- switch(ncpStat, mode=cmpnt$ncp_mode, median=cmpnt$ncp_median,mean=cmpnt$ncp,pct90=cmpnt$ncp_pct90,max=sum(!is.nan(cp)))
ncp <<- base::round(ncp)
ncpPr <<- cmpnt$ncpPr
cpPr <<- cmpnt$cpPr
cpChange<<- cmpnt$cpAbruptChange
HasChangePoint <<- !is.null(x$trend$cp);
Prob <<- cmpnt$cpOccPr;
Prob2Prob <<- c(Prob,Prob-Prob)
}
###########################################################
Amp = 0
AmpSD = 0
get.S = function(){
Y <<- list()
SD <<- list()
CI <<- list()
Amp <<- list()
AmpSD <<- list()
for (i in 1:q){
Y[[i]] <<-x$season$Y[[i]];
SD[[i]] <<-c(Y[[i]] -x$season$SD[[i]], rev(Y[[i]] +x$season$SD[[i]]));
if ( !is.null(x$season$CI) ) CI[[i]] <<-c(x$season$CI[[i]][,1], rev(x$season$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
if (hasAmp){
Amp[[i]] <<-x$season$amp[[i]]
AmpSD[[i]] <<-c(Amp[[i]]-x$season$ampSD[[i]], rev(Amp[[i]]+x$season$ampSD[[i]]));
}
}
Order <<-x$season$order;
HasChangePoint <<- !is.null(x$season$cp);
}
get.scp = function(){
cmpnt = x$season
cp <<-cmpnt$cp;
cpCI <<-cmpnt$cpCI;
ncp <<-switch(ncpStat, mode=cmpnt$ncp_mode, median=cmpnt$ncp_median,mean=cmpnt$ncp,pct90=cmpnt$ncp_pct90,max=sum(!is.nan(cp)))
ncp <<-base::round(ncp)
ncpPr <<-cmpnt$ncpPr
cpPr <<-cmpnt$cpPr
cpChange<<-cmpnt$cpAbruptChange
Prob <<-x$season$cpOccPr;
Prob2Prob <<-c(Prob,Prob-Prob)
HasChangePoint <<- !is.null(x$season$cp);
}
get.O = function(){
Y <<-list()
SD <<-list()
for (i in 1:q){
Y[[i]] <<-x$outlier$Y[[i]] ;
SD[[i]] <<-c(Y[[i]] -x$outlier$SD[[i]] , rev(Y[[i]] +x$outlier$SD[[i]] ));
if ( !is.null(x$outlier$CI) ) CI[[i]] <<-c(x$outlier$CI[[i]][,1], rev(x$outlier$CI[[i]][,2]))
else CI[[i]] <<-SD[[i]]
}
}
get.ocp = function(){
#cp <<-x$season$cp;
#cpCI <<-x$season$cpCI;
#ncp <<-sum(!is.nan(cp))
#ncpPr<<-x$season$ncpPr
#cpPr <<-x$season$cpPr
Prob <<-x$outlier$cpOccPr;
Prob2Prob <<-c(Prob,Prob-Prob)
}
YSIDE = 2;
padj = -1.2
tcl = 0.2
plotId = 1;
plot.st=function(col=c(0.2,0.2,0.2), ylabel){
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
if (hasData){
Yall=c(YtsSD[[i]],data[[i]])
plot( c(t[1],t[N]), c(min(Yall,na.rm = TRUE),max(Yall,na.rm = TRUE)), type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
} else {
plot( t2t, YtsSD[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
}
### polygon(t2t, YtsSD[[i]], col = rgb(.5,0.5,0.5,.5), border = NA);
if (hasData){
points( t, data[[i]], type = 'p', col='#777777');
}
points( t, Yts[[i]], type = 'l', lwd=2,col=rgb(col[1],col[2],col[3]));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.y =function(col=c(0,1,0), ylabel){
alpha=0.1
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
if (hasData && !hasSeason){
Yall=c(CI[[i]],data[[i]])
plot( c(t[1],t[N]), c(min(Yall,na.rm = TRUE),max(Yall,na.rm = TRUE)), type = 'n',ann=FALSE, xaxt='n', yaxt='n');
} else {
plot( t2t, CI[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
}
if (hasData && !hasSeason){
points(t, data[[i]], type = 'p', col='#777777');
}
### polygon(t2t, CI[[i]], col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Y[[i]], type = 'l',col=rgb(col[1],col[2],col[3] ) );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.prob <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
alpha=0.2
plot( c(t2t[1],t2t), c(0.22,Prob2Prob),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
polygon(t2t, Prob2Prob, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Prob, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.order <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
maxOrder=max(1.01,max(Order));
plot( c(t[1],t[length(t)]), c(0,maxOrder),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
points( t, Order, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
yext =par('usr')
if (HasChangePoint) {
for (i in 1 : ncp) {
points(c(cp[i], cp[i]), yext[3:4],type = 'l', lty = 'dashed', col = 'grey10',lwd =1, bg = 'grey10');
}
}
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.amp =function(col=c(0,1,0), ylabel){
alpha=0.5
plot( t, Amp, type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, AmpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Amp, type = 'l',col=rgb(col[1],col[2],col[3] ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.slp =function(col=c(0,1,0), ylabel){
alpha=0.5
plot( t2t, SlpSD, type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
polygon(t2t, SlpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
points( t, Slp, type = 'l',col=rgb(0,0,0,0.8 ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.slpsgn = function(col=c(0,1,0), ylabel){
alpha=0.5
SlpSignNeg = 1-SlpSignPos-SlpSignZero
plot( c(t[1],t[length(t)]), c(0,1), type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, SlpSD, col = rgb(col[1],col[2],col[3],alpha), border = NA);
y2y = c(t-t, rev(SlpSignNeg)); polygon(t2t, y2y, col=rgb(0,0,1,alpha),border=NA)
y2y = c(SlpSignNeg, rev(1-SlpSignPos)); polygon(t2t, y2y, col=rgb(0,1,0,alpha),border=NA)
y2y = c(1-SlpSignPos, t-t+1); polygon(t2t, y2y, col=rgb(1,0,0,alpha),border=NA)
points( t, t-t+0.5 , type = 'l', lty=2, col=rgb(col[1],col[2],col[3] ));
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.o =function(col=c(0,1,0), ylabel){
alpha=0.5
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
plot( t2t, CI[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
#polygon(t2t, CI, col = rgb(col[1],col[2],col[3],alpha), border = NA);
#points( t, Y, type = 'l',col='#333333');
segments(t,t-t,t,Y[[i]], col= rgb(col[1],col[2],col[3]))
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.oprob <- function( col=c(0,1,0), ylabel ){
subplot(plotId,heights )
plotId <<- plotId+1;
alpha=0.2
plot( c(t2t[1],t2t), c(0.22,Prob2Prob),type = 'n', ann=FALSE, xaxt='n', yaxt='n');
# polygon(t2t, Prob2Prob, col = rgb(col[1],col[2],col[3],alpha), border = NA);
# points( t, Prob, col = rgb(col[1],col[2],col[3]) , lwd = 1,type = 'l' );
segments(t,t-t,t,Prob, col= rgb(col[1],col[2],col[3]),lwd=1.5)
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
plot.error=function(col=c(0.2,0.2,0.2), ylabel){
for (i in 1:q) {
subplot(plotId,heights )
plotId <<- plotId+1;
plot( t, Yerr[[i]], type = 'n',ann = FALSE, xaxt = 'n', yaxt = 'n');
lines(t, t-t, col= rgb(col[1],col[2],col[3]))
segments(t,t-t,t,Yerr[[i]], col= rgb(col[1],col[2],col[3]))
axis(1, labels = FALSE, padj = padj, tcl = tcl);
axis(YSIDE, labels = TRUE, padj = ifelse(YSIDE==2,-padj,padj), tcl = tcl)
YSIDE<<-(YSIDE==2)*4+(YSIDE==4)*2
title(ylab = ylab[plotId-1],cex=cex.lab, mgp = c(1.25, 1.25, 0));
}
}
plot.new()
opar= par(c('oma','mar','fig'))
par(oma=c(0,0,0,0))
par(mar=c(0,2,0,2))
par(fig=c(0,1,0,1))
lm= strheight('a',units='figure')*dev.size()[2]/dev.size()[1]*4.5
rm= strheight('a',units='figure')*dev.size()[2]/dev.size()[1]*2.3
bm= strheight('a',units='figure')*4.0
tm= strheight('a',units='figure')*4.0
#######################################################
# Create a subplot given the relative heights of the vertical plots
########################################################
subplot <- function(i, heights){
Hmiddle =1-tm-bm
h = (heights/sum(heights))*Hmiddle
hcum = cumsum(h)
plt = c(lm, 1-rm, 1-tm-hcum[i], 1-tm-hcum[i]+h[i])
par(plt=plt,new=TRUE)
}
for (i in 1:length(vars)){
ytitle = ylab[i]
var = tolower(vars[i])
clr = as.numeric(col2rgb( col[i] ))/255
if (var=='st') { get.Yts(); plot.st(clr, ytitle) }
if (var=='s' ) { get.S(); get.scp(); plot.y(clr, ytitle) }
if (var=='t' ) { get.T(); get.tcp(); plot.y(clr, ytitle) }
if (var=='scp') { get.S(); get.scp(); plot.prob(clr, ytitle) }
if (var=='tcp') { get.T(); get.tcp(); plot.prob(clr, ytitle) }
if (var=='sorder'){ get.S(); get.scp(); plot.order(clr, ytitle) }
if (var=='torder'){ get.T(); get.tcp(); plot.order(clr, ytitle) }
if (var=='samp') { get.S(); plot.amp(clr, ytitle) }
if (var=='tslp') { get.T(); plot.slp(clr, ytitle) }
if (var=='slpsgn'){ get.T(); plot.slpsgn(clr, ytitle) }
if (var=='o') { get.O(); plot.o(clr, ytitle) }
if (var=='ocp') { get.ocp(); plot.oprob(clr, ytitle) }
if (var=='error') { get.Yts(); plot.error(clr, ytitle) }
#title(ylab = ytitle,cex=cex.lab, mgp = c(1.25, 1.25, 0));
if(i==1){
#mtext(main,line=.5,cex=cex.main,font=2)
}
if(i==nPlots){
title(xlab = xlab ,mgp = c(1.2, 1.2, 0));
if (sum(class(t)=='Date') >0 ){
axis.Date(1, t, labels = TRUE, padj =padj, tcl =tcl,cex.axis=1);
} else{
axis( 1, labels = TRUE, padj =padj, tcl =tcl,cex.axis=1);
}
}
}
par(opar)
#text(mean(time) / 2.0, 'NO SEASONAL COMPONET.IGNORE THIS POLOT');
}
#plot.mrbeast(o,vars=c('st','s','scp','t','tcp','o','ocp'))
Try the Rbeast package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.