beast.irreg  R Documentation 
A Bayesian model averaging algorithm called BEAST to decompose time series or 1D sequential data into individual components, such as abrupt changes, trends, and periodic/seasonal variations. BEAST is useful for changepoint detection (e.g., breakpoints or structural breaks), nonlinear trend analysis, time series decomposition, and time series segmentation.
beast.irreg(
y,
time,
deltat = NULL,
season = c("harmonic", "svd", "dummy", "none"),
period = NULL,
scp.minmax = c(0,10), sorder.minmax = c(0,5),
tcp.minmax = c(0,10), torder.minmax = c(0,1),
sseg.min = NULL, sseg.leftmargin = NULL, sseg.rightmargin = NULL,
tseg.min = NULL, tseg.leftmargin = NULL, tseg.rightmargin = NULL,
method = c('bayes', 'bic', 'aic', 'aicc', 'hic',
'bic0.25', 'bic0.5', 'bic1.5', 'bic2' ),
detrend = FALSE,
deseasonalize = FALSE,
mcmc.seed = 0,
mcmc.burnin = 200,
mcmc.chains = 3,
mcmc.thin = 5,
mcmc.samples = 8000,
precValue = 1.5,
precPriorType = c('componentwise', 'uniform', 'constant', 'orderwise'),
hasOutlier = FALSE,
ocp.minmax = c(0,10),
print.param = TRUE,
print.progress = TRUE,
print.warning = TRUE,
quiet = FALSE,
dump.ci = FALSE,
dump.mcmc = FALSE,
gui = FALSE,
... )
y 
a vector for an irregular or unordered time series. Missing values such as NA and NaN are allowed.
If 
time 
a vector of the same length as

deltat 
a number or a string to specify a time interval for aggregating the irregular 
season 
characters (default to 'harmonic'); specify if

period 
numeric or string. Specify the period for the periodic/seasonal component in 
scp.minmax 
a vector of 2 integers (>=0); the min and max number of seasonal changepoints (scp) allowed in segmenting the seasonal component. 
sorder.minmax 
a vector of 2 integers (>=1); the min and max harmonic orders considered to fit the seasonal component. 
tcp.minmax 
a vector of 2 integers (>=0); the min and max number of trend changepoints (tcp) allowed in segmenting the trend component. If the min and max changepoint numbers are equal, BEAST assumes a constant number of changepoints and won't infer the posterior probability of the number of changepoints for the trend, but it still estimates the occurrence probability of the changepoints over time (i.e., the most likely times at which these changepoints occur in the trend). If both the min and max numbers are set to 0, no changepoints are allowed; then a global polynomial trend is used to fit the trend component, but still, the most likely polynomial order will be inferred if torder.minmax[1] is not equal to torder.minmax[2]. 
torder.minmax 
a vector of 2 integers (>=0); the min and max orders of the polynomials considered to fit the trend component. The 0th order corresponds to a constant term/a flat line and the 1st order is a line. If 
sseg.min 
an integer (>0); the min segment length allowed between two neighboring season changepoints. That is, when fitting a piecewise harmonic seasonal model, two changepoints are not allowed to occur within a time window of length 
sseg.leftmargin 
an integer (>=0); the number of leftmost data points excluded for seasonal changepoint detection. That is, when fitting a piecewise harmonic seasonal model, no changepoints are allowed in the starting window/segment of length 
sseg.rightmargin 
an integer (>=0); the number of rightmost data points excluded for seasonal changepoint detection. That is, when fitting a piecewise harmonic seasonal model, no changepoints are allowed in the ending window/segment of length 
tseg.min 
an integer (>0); the min segment length allowed between two neighboring trend changepoints. That is, when fitting a piecewise polynomial trend model, two changepoints are not allowed to occur within a time window of length 
tseg.leftmargin 
an integer (>=0); the number of leftmost data points excluded for trend changepoint detection. That is, when fitting a piecewise polynomial trend model, no changepoints are allowed in the starting window/segment of length 
tseg.rightmargin 
an integer (>=0); the number of rightmost data points excluded for trend changepoint detection. That is, when fitting a piecewise polynomial trend model, no changepoints are allowed in the ending window/segment of length 
method 
a string (default to 'bayes'); specify the method for formulating model posterior probability.

detrend 
logical; If 
deseasonalize 
logical; If 
mcmc.seed 
integer (>=0); the seed for the random number generator used for Monte Carlo Markov Chain (mcmc). If 
mcmc.chains 
integer (>0); the number of MCMC chains. 
mcmc.thin 
integer (>0); a factor to thin chains (e.g., if thinningFactor=5, samples will be taken every 3 iterations) 
mcmc.burnin 
integer (>0); the number of burnin samples discarded at the start of each chain 
mcmc.samples 
integer (>=0); the number of samples collected per MCMC chain. The total number of iterations is 
precValue 
numeric (>0); the hyperparameter of the precision prior; the default value is 1.5. 
precPriorType 
characters. It takes one of 'constant', 'uniform', 'componentwise' (the default), and 'orderwise'. Below are the differences between them.

hasOutlier 
boolean; if true, fit a model with an outlier component that refers to potential spikes or dips at isolated data points: Y = trend + outlier + error if season='none',and Y = trend + season + outlier + error if season ~= 'none'. 
ocp.minmax 
a vector of 2 integers (>=0); the min and max numbers of outliertype changepoints (ocp) allowed in the time seriestrend component. Ocp refers to spikes or dips at isolated times that can't be modeled as trends or seasonal terms. 
print.param 
boolean. If 
print.progress 
boolean;If 
print.warning 
boolean;If 
quiet 
boolean. If 
dump.ci 
boolean; If 
dump.mcmc 
boolean; If 
gui 
boolean. If 
... 
additional parameters. There are many more settings for the implementation but not made available in the beast() interface; please use the function 
The output is an object of class "beast". It is a list, consisting of the following variables. In the explanations below, we assume the input y
is a single time series of length N
:
time 
a vector of size 
data 
a vector, matrix, or 3D array; this is a copy of the input 
marg_lik 
numeric; the average of the model marginal likelihood; the larger marg_lik, the better the fitting for a given time series. 
R2 
numeric; the Rsquare of the model fitting. 
RMSE 
numeric; the RMSE of the model fitting. 
sig2 
numeric; the estimated variance of the model error. 
trend 
a list object consisting of various outputs related to the estimated trend component:

season 
a list object consisting of various outputs related to the estimated seasonal/periodic component:

The three functions beast
(), beast.irreg
(), and beast123
() are essentially the same BEAST algorithm but with different APIs. There is a onetoone correspondence between the parameters for beast() and beast.irreg() and the 'metadata', 'prior','mcmc', and 'extra' objects in the beast123() interface. Examples are:
start <> metadata$startTime
deltat <> metadata$deltaTime
deseasonalize <> metadata$deseasonalize
hasOutlier <> metadata$hasOutlierCmpnt
scp.minmax[1] <> prior$seasonMinOrder
scp.minmax[2] <> prior$seasonMaxOrder
sseg.min <> prior$seasonMinSepDist
tcp.torder[1] <> prior$trendMinOrder
tseg.leftmargin <> prior$trendLeftMargin
mcmc.seed <> mcmc$seed
dump.ci <> extra$computeCredible
Experts should use the the beast123 function.
Zhao, K., Wulder, M.A., Hu, T., Bright, R., Wu, Q., Qin, H., Li, Y., Toman, E., Mallick, B., Zhang, X. and Brown, M., 2019. Detecting changepoint, trend, and seasonality in satellite time series data to track abrupt changes and nonlinear dynamics: A Bayesian ensemble algorithm. Remote Sensing of Environment, 232, p.111181 (the beast algorithm paper).
Zhao, K., Valle, D., Popescu, S., Zhang, X. and Mallick, B., 2013. Hyperspectral remote sensing of plant biochemistry using Bayesian model averaging with variable and band selection. Remote Sensing of Environment, 132, pp.102119 (the Bayesian MCMC scheme used in beast).
Hu, T., Toman, E.M., Chen, G., Shao, G., Zhou, Y., Li, Y., Zhao, K. and Feng, Y., 2021. Mapping finescale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250261(a beast application paper).
beast
, beast123
, minesweeper
, tetris
, geeLandsat
library(Rbeast)
######################################################################################
# Note that the BEAST algorithm is currently implemented to handle only regular time
# series. 'beast.irreg' accepts irregular time series but internally it aggregates them
# into regular ones prior to applying the BEAST model. For the aggregation, both the
# "time" and "deltat" args are needed to specify individual times of data points and the
# regular time interval desired. If there is a cyclic componet, 'period' should also be given;
# if not, a possible value is guessed via autocorrelation
######################################################################################
# 'ohio' is a data.frame on an irregular Landsat time series of reflectances & ndvi
# (e.g., surface greenness) at an Ohio site. It has multiple columns of alternative date
# formats, such as year, month, day, doy (date of year), rdate (R's date class), and
# time (fractional year)
data(ohio)
str(ohio)
plot(ohio$rdate, ohio$ndvi,type='o') # ndvi is irregularly spaced and unordered in time
######################################################################################
# Below, 'time' is given as numeric values, which can be of any arbitray unit. Although
# here 1/12 can be interepreted as 1/12 year or 1 month, BEAST itself doesn't care about
# the time unit. So, the unit of 1/12 is irrelevant for BEAST. 'freq' or 'period' is missing
# and a guess of it is used.
o=beast.irreg(ohio$ndvi, time=ohio$time,deltat=1/12)
plot(o)
print(o)
######################################################################################
# Aggregrate the time series at a monthly interval (deltat=1/12) and explictly provide
# the 'freq' or 'period' arg
o=beast.irreg(ohio$ndvi, time=ohio$time,deltat=1/12, period=1.0)
#o=beast.irreg(ohio$ndvi, time=ohio$time,deltat=1/12, freq =12)
## Not run:
######################################################################################
# Aggregate the time series at a halfmonthly time interval, and the 'freq' becomes 24
# while the period is still 1. That is, PERIOD (1.0)=deltat(1/24) X freq (24)
o=beast.irreg(ohio$ndvi, time=ohio$time,deltat=1/24, freq = 24)
#o=beast.irreg(ohio$ndvi, time=ohio$time,deltat=1/24, period = 1)
######################################################################################
# 'time' is given as R's dates. The unit is YEAR. 1/12 refers to 1/12 year or 1 month
o=beast.irreg(ohio$ndvi, time=ohio$rdate,deltat=1/12)
######################################################################################
# 'time' is given as data strings. The unit is YEAR. 1/12 refers to 1/12 year or 1 month
o=beast.irreg(ohio$ndvi, time=ohio$datestr1,deltat=1/12) #"LT419840327" (YYYYMMDD)
o=beast.irreg(ohio$ndvi, time=ohio$datestr2,deltat=1/12) #"LT41984087ndvi" (YYYYDOY)
o=beast.irreg(ohio$ndvi, time=ohio$datestr3,deltat=1/12) #"1984,, 3/ 27" (YYYY M D)
######################################################################################
# 'time' is given as data strings, with a format specifier
TIME =list()
TIME$datestr = ohio$datestr1
TIME$strfmt = "LT4YYYYMMDD" # "LT419840327"
o=beast.irreg(ohio$ndvi, time=TIME,deltat=1/12)
TIME =list()
TIME$datestr = ohio$datestr2
TIME$strfmt = "LT4YYYYDOYndvi" # LT41984087ndvi
o=beast.irreg(ohio$ndvi, time=TIME,deltat=1/12)
######################################################################################
# 'time' is given as a list object
TIME = list()
TIME$year = ohio$Y
TIME$month = ohio$M
TIME$day = ohio$D
o=beast.irreg(ohio$ndvi, time=TIME,deltat=1/12)
TIME = list()
TIME$year = ohio$Y
TIME$doy = ohio$doy
o=beast.irreg(ohio$ndvi, time=TIME, deltat=1/12)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.