beast.irreg: Bayesian decomposition of irregular time series for...

In Rbeast: Bayesian Change-Point Detection and Time Series Decomposition

Bayesian decomposition of irregular time series for changepoints, trend, and seasonality

Description

beast.irreg applies the BEAST (Bayesian Estimator of Abrupt change, Seasonal change, and Trend) model to irregular or unordered time series or 1D sequential data. BEAST is a Bayesian model averaging algorithm for decomposing time series or other 1D sequential data into individual components, including abrupt changes, trends, and periodic/seasonal variations. It is useful for changepoint detection (e.g., breakpoints or structural breaks), nonlinear trend analysis, time series decomposition, and time series segmentation. beast123 is a low-level interface to BEAST.

Internally, irregular observations are first aggregated into an evenly spaced (regular) time series at a user-specified time interval, and then decomposed into individual components such as abrupt changes, nonlinear trends, and periodic/seasonal variations.

Usage

  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, 
       s.complexfct   = 0.0,
       t.complexfct   = 0.0,		   
       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,
       ...
  )

Arguments

y	y is a numeric vector for an irregular or unordered time series. Missing values such as NA and NaN are allowed. If y is evenly spaced in time (regular), use beast instead. If y is a matrix or 3D array (e.g., a stack of images) consisting of multiple regular or irregular time series, use beast123 instead. If y is a list of multiple time series, the multivariate version of the BEAST algorithm is invoked to jointly decompose the series and detect common changepoints. This feature is currently experimental and under development. See ohio for a working example.
time	time is a vector (or list) specifying the times of each observation in y. Its length must match the time dimension of y. It can be numeric, Date, character, or a list of date components. Supported formats include: Numeric vector For example c(1984.23, 1984.27, 1984.36, ...). The unit of time is irrelevant to BEAST as long as it is used consistently with deltat and period. R Date vector For example as.Date(c("1984-03-27", "1984-04-10", "1984-05-12", ...)). Character vector of date strings For example: c("1984-03-27", "1984-04-10", "1984-05-12") c("1984/03/27", "1984,04,10", "1984 05 12") (delimiters may differ as long as the year-month-day order is consistent) c("LT4-1984-03-27", "LT4-1984-04-10", "LT4-1984+05,12") c("LT4-1984087ndvi", "LT4-1984101ndvi", "LT4-1984133ndvi") c("1984,,abc 3/ 27", "1984,,ddxfdd 4/ 10", "ggd1984,, 5/ ttt 12") BEAST uses several heuristics to parse date strings without an explicit format specifier and may fail for ambiguous patterns (e.g., in "LC8-2020-09-20-1984" it is unclear whether 2020 or 1984 is the year). For robust parsing, consider using a list with an explicit format (see the next). List with date strings and format A list of the form time = list(datestr = ..., strfmt = "..."), where time$datestr is a character vector of date strings and time$strfmt is a format specifier describing how to extract year, month, day, or day-of-year. Three classes of formats are supported: (a) Fixed positions for year/month/day. For example, to extract 2001/12/02 etc. from time$datestr = c("P23R34-2001.1202333xd", "O93X94-2002.1108133fd", "TP3R34-2009.0122333td") use time$strfmt = "P23R34-yyyy.mmdd333xd", where yyyy, mm, and dd mark the positions of year, month, and day; all other characters are wildcards. (b) Fixed positions for year and day-of-year (doy). For example, to extract years and days-of-year from strings like "P23R342001888045" use time$strfmt = "123123yyyy888doy", where yyyy and doy mark year and day-of-year; all other positions are wildcards. doy must be three digits. (c) Patterns based on separators between Y/M/D. For example, to extract 2002/12/02 from "2002,12/02", " 2002 , 12/2", "2002,12 /02 " use time$strfmt = "Y,M/D"; whitespace is ignored. To extract 2002/12/02 from "2--12, 2012", use time$strfmt = "D--M,Y". List of date components A list specifying dates component-wise, using either time$year, time$month, and time$day, or time$year and time$doy where each element is a vector of the same length as the time dimension of y. For the doy representation, values must lie between 1 and 365/366.
deltat	Numeric or character; the time interval used to aggregate the irregular time series into a regular one. The BEAST model is formulated for regular (evenly spaced) data; therefore beast.irreg first re-bins the irregular time series into regular intervals of length deltat. If deltat is missing, a heuristic guess is used. The unit of deltat must be consistent with time: If time is numeric, the unit of deltat is arbitrary but must be consistent with time. If time is a Date vector or date strings, the default unit of deltat is fractional years. To specify the unit explicitly, supply a character string, for example "7 days", "7d", "1/2 months", "1 mn", "1.0 year", or "1y".
season	Character (default "harmonic"); indicates whether y has a periodic component, and how it should be modeled. Allowed values are: "none": y is trend-only; no periodic component is modeled. Seasonal arguments (e.g., sorder.minmax, scp.minmax, sseg.min) are ignored. "harmonic": y has a periodic/seasonal component. Here, ‘season’ is used broadly to mean any periodic variation. The period is a known constant supplied via period and is not estimated as a model parameter. The seasonal component is modeled as a harmonic curve (sum of sines and cosines). "dummy": As in "harmonic", but the periodic component is modeled via non-parametric dummy bases. The harmonic-order argument sorder.minmax is irrelevant and ignored. "svd": (experimental) As in "harmonic", but the seasonal component is modeled as a linear combination of basis functions derived via singular value decomposition (SVD). These SVD-based bases can be more parsimonious than the standard harmonic bases and can improve sensitivity to subtle changepoints.
period	Numeric or character. The period of the seasonal/periodic component in y. Required only when a periodic component is present (i.e., season is "harmonic", "svd", or "dummy"), and ignored if season = "none". The period must have units consistent with deltat, and satisfies period = deltat * freq, where freq is the number of samples per period. The historical freq argument (number of data points per period) is still supported via ..., but is deprecated; if both period and freq are supplied, period takes precedence. period (or equivalently freq) is not itself a BEAST model parameter and must usually be specified by the user. If period is missing, BEAST attempts to guess it via autocorrelation before fitting the model. If period <= 0, season = "none" is assumed and a trend-only model is fitted. To specify units explicitly for date-based time axes, use a string such as "1.0 year", "12 months", "365d", or "366 days".
scp.minmax	Integer vector of length 2 (>= 0); minimum and maximum numbers of seasonal changepoints (SCPs) allowed in the seasonal component. Used only when a seasonal component is present (i.e., season != "none") and ignored otherwise. If scp.minmax[1] == scp.minmax[2], BEAST assumes a fixed number of seasonal changepoints and does not infer the posterior on the number of changepoints, but still estimates when those changepoints occur. If both values are 0, no seasonal changepoints are allowed; a global harmonic (or SVD/dummy) model is used, but the most likely harmonic order is still inferred if sorder.minmax[1] != sorder.minmax[2].
sorder.minmax	Integer vector of length 2 (>= 1); minimum and maximum harmonic orders considered for the seasonal component. Used only if season = "harmonic" or "svd" and ignored for trend-only data or when season = "dummy". If sorder.minmax[1] == sorder.minmax[2], BEAST assumes a fixed harmonic order and does not infer the posterior distribution of harmonic orders.
tcp.minmax	Integer vector of length 2 (>= 0); minimum and maximum numbers of trend changepoints (TCPs) allowed in the trend component. If tcp.minmax[1] == tcp.minmax[2], BEAST assumes a fixed number of trend changepoints and does not infer the posterior on the number of trend changepoints, but still estimates their occurrence probabilities over time. If both values are 0, no trend changepoints are allowed; a global polynomial trend is used, but the most likely polynomial order is still inferred if torder.minmax[1] != torder.minmax[2].
torder.minmax	Integer vector of length 2 (>= 0); minimum and maximum polynomial orders considered for the trend component. Order 0 corresponds to a constant (flat) trend; order 1 is a line. If torder.minmax[1] == torder.minmax[2], BEAST assumes a fixed polynomial order and does not infer the posterior distribution of trend orders.
sseg.min	Integer (> 0); minimum length (in data points) of a segment between two neighboring seasonal changepoints. When fitting a piecewise seasonal model, no two seasonal changepoints are allowed to occur within sseg.min time steps. sseg.min is unitless (number of time intervals); in the original time unit the window length is sseg.min * deltat. If NULL, a default based on the implied frequency is used.
sseg.leftmargin	Integer (>= 0); number of leftmost data points excluded from seasonal changepoint detection. No seasonal changepoints are allowed in the initial window of length sseg.leftmargin. sseg.leftmargin is unitless (number of samples). In the original time unit, the excluded window length is sseg.leftmargin * deltat. If NULL, defaults to sseg.min.
sseg.rightmargin	Integer (>= 0); number of rightmost data points excluded from seasonal changepoint detection. No seasonal changepoints are allowed in the ending window of length sseg.rightmargin. sseg.rightmargin is unitless, and the corresponding time window is sseg.rightmargin * deltat. If NULL, defaults to sseg.min.
tseg.min	Integer (> 0); minimum length (in data points) of a segment between two neighboring trend changepoints. When fitting a piecewise polynomial trend, no two trend changepoints are allowed within a window of length tseg.min. tseg.min is unitless; in the original time unit the window is tseg.min * deltat. If NULL, a default based on the implied frequency is used when a cyclic component is present.
tseg.leftmargin	Integer (>= 0); number of leftmost data points excluded from trend changepoint detection. No trend changepoints are allowed in the starting window of length tseg.leftmargin. tseg.leftmargin is unitless; the excluded time window is tseg.leftmargin * deltat. If NULL, defaults to tseg.min.
tseg.rightmargin	Integer (>= 0); number of rightmost data points excluded from trend changepoint detection. No trend changepoints are allowed in the ending window of length tseg.rightmargin. tseg.rightmargin is unitless; the excluded time window is tseg.rightmargin * deltat. If NULL, defaults to tseg.min.
s.complexfct	Numeric (default 0.0). A hyperprior parameter–newly added in Version 1.02– controlling the complexity of the seasonal curve (i.e., the number of seasonal changepoints). A prior of the form p(k) \propto \exp[\lambda (k+1)] is placed on the number of seasonal changepoints k, where \lambda is seasonComplexityFactor. Setting \lambda = 0 yields a non-informative prior p(k) \propto 1.0 where all model dimensions are equally likely a priori. Users may tune seasonComplexityFactor in the range [-20, 20] or an even wider range: negative values (e.g., \lambda = -15.9) favor fewer changepoints (simpler seasonal curves), whereas positive values (e.g., \lambda = 5.76) favor more changepoints (more complex curves).
t.complexfct	Numeric scalar (default 0.0); analogous to s.complexfct but for the trend component and the number of trend changepoints.
method	Character (default "bayes"); specifies how model posterior probabilities are formulated or approximated: "bayes": full Bayesian formulation as in Zhao et al. (2019). "bic": BIC approximation, bic = n * ln(SSE) + k * ln(n), where k and n are the numbers of parameters and data points. "aic": AIC approximation, aic = n * ln(SSE) + 2k. "aicc": corrected AIC, aicc = aic + (2k^2 + 2k) / (n - k - 1). "hic": Hannan–Quinn information criterion, hic = n * ln(SSE) + 2k * ln(ln(n)). "bic0.25": BIC-type approximation adopted from Kim et al. (2016) <doi:10.1016/j.jspi.2015.09.008>, bic0.25 = n * ln(SSE) + 0.25 k * ln(n), with a weaker penalty on model complexity. "bic0.50": same as above but with penalty factor 0.50. "bic1.5": same as above but with penalty factor 1.5. "bic2": same as above but with penalty factor 2.0.
detrend	Logical; if TRUE, a global trend is first fitted and removed from the series before running BEAST, and then added back to the final result. This can help when the global trend dominates other components.
deseasonalize	Logical; if TRUE, a global seasonal model is first fitted and removed from the series before running BEAST, then added back to the final result. Ignored when season = "none".
mcmc.seed	Integer (>= 0); seed for the random number generator used in the MCMC sampler. If mcmc.seed = 0, an arbitrary seed is chosen and results vary across runs. Using a fixed non-zero seed improves reproducibility on the same machine; results may still differ across architectures because the underlying random number generator can depend on CPU instruction sets.
mcmc.chains	Integer (> 0); number of parallel MCMC chains.
mcmc.thin	Integer (> 0); thinning factor for the chains (e.g., if mcmc.thin = 5, every 5th iteration is retained).
mcmc.burnin	Integer (> 0); number of burn-in samples discarded at the beginning of each chain.
mcmc.samples	Integer (>= 0); number of post-burn-in samples collected per MCMC chain. The total number of iterations is (mcmc.burnin + mcmc.samples * mcmc.thin) * mcmc.chains.
precValue	Numeric (> 0); hyperparameter for the precision prior. It is directly used only when precPriorType = "constant"; otherwise it serves as an initial value that is updated within MCMC (see below).
precPriorType	Character; one of "constant", "uniform", "componentwise" (default), or "orderwise". These options control how precision parameters for model coefficients are treated: "constant": a single precision parameter is fixed to precValue. The result may be sensitive to precValue. "uniform": a single precision parameter is treated as a random variable; its initial value is precValue, but its posterior is inferred via MCMC, making results less sensitive to the initial choice. "componentwise": separate precision parameters are used for different components (e.g., one for season, another for trend), initialized by precValue and updated via MCMC. "orderwise": separate precision parameters are used not only for components but also for individual orders within each component, again initialized by precValue and inferred via MCMC.
hasOutlier	Logical; if TRUE, fits a model with an additional outlier component capturing spikes or dips at isolated data points: Y = trend + outlier + error when season = "none", Y = trend + season + outlier + error otherwise. Outliers are modeled as changepoints that cannot be captured by trend or seasonal terms.
ocp.minmax	Integer vector of length 2 (>= 0); minimum and maximum numbers of outlier-type changepoints (OCPs) allowed in the series. OCPs correspond to isolated spikes/dips that are not captured by trend or seasonality.
print.param	Logical; if TRUE, the full list of internal BEAST parameters (organized as metadata, prior, mcmc, and extra) is printed before MCMC. The internal naming differs slightly from the beast.irreg interface, but there is a one-to-one mapping (e.g., prior$trendMinSepDist <- tseg.min). See beast123 or View(beast) for details.
print.progress	Logical; if TRUE, prints a progress bar during model fitting.
print.warning	Logical; if TRUE, prints warning messages.
quiet	Logical; if TRUE, suppresses all console output (overrides other printing options).
dump.ci	Logical; if TRUE, computes explicit credible intervals (i.e., out$season$CI, out$trend$CI) for the estimated seasonal and trend components. This requires sorting and can be time-consuming. If only symmetric intervals are needed, use the standard deviations (out$trend$SD, out$season$SD) and set dump.ci = FALSE.
dump.mcmc	Logical; if TRUE, dumps individual MCMC samples for further custom analysis.
gui	Logical; if TRUE, launches a GUI window that animates MCMC sampling in the model space (numbers and timings of seasonal and trend changepoints). This experimental GUI is currently available only on 64-bit Windows, not on 32-bit Windows, macOS, or Linux.
...	Additional implementation-level arguments passed to the underlying beast123() engine. For full control over advanced settings, use beast123 directly.

Value

An object of class "beast", i.e., a list with components similar to those returned by beast and beast123. Below we assume the input y is a single time series of length N:

time	Numeric vector of length N; the regularized time axis after aggregating the original irregular series. If no explicit time is available, it defaults to 1:N.
data	A vector, matrix, or 3D array containing a copy of the regularized input data if extra$dumpInputData = TRUE. Otherwise NULL. If the original input was irregular, this field stores the aggregated regular series at the time interval specified by metadata$deltaTime.
marg_lik	Numeric; the (average) model marginal likelihood. Larger values correspond to better fits for a given time series (e.g., -1 is better than -10; 10 is better than -1).
R2	Numeric; coefficient of determination (R^2) for the fitted model.
RMSE	Numeric; root mean squared error of the fitted model.
sig2	Numeric; estimated variance of the model residuals.
trend	List containing summaries for the estimated trend component: ncp: Mean number of trend changepoints across sampled models. Because individual models have varying numbers of changepoints, additional summaries are provided: ncp_mode, ncp_median, ncp_pct90. For example, if mcmc$samples = 10 yields changepoint counts c(0, 2, 4, 1, 1, 2, 7, 6, 6, 1), then the mean is 3.1, median 2.5, mode 1, and 90th percentile 6.5. ncp_mode: Mode of the posterior distribution of the number of trend changepoints. ncp_median: Median of the posterior distribution of the number of trend changepoints. ncp_pct90: 90th percentile of the posterior distribution of the number of trend changepoints. ncpPr: Vector of length tcp.max + 1 where tcp.max = tcp.minmax[2]; gives the probability of having 0, 1, ..., tcp.max trend changepoints. For example, ncpPr[1] is the probability of having no trend changepoint; ncpPr[i] is the probability of having i - 1 changepoints. cpOccPr: Numeric vector of length N; the marginal probability at each time index of being a trend changepoint. Plotting cpOccPr yields a continuous curve of "probability of being a changepoint". Note: a higher peak in this curve reflects a higher probability at that specific time index, but not necessarily a higher probability of a changepoint in a neighborhood around that time. For instance, a window c(0, 0, 0.5, 0, 0) (peak 0.5, sum 0.5) can be less likely to contain a changepoint than c(0.1, 0.2, 0.21, 0.2, 0.1) (peak 0.21, sum 0.71). order: Numeric vector of length N; average polynomial order of the trend at each time index, averaged over sampled piecewise polynomial trends. It is not necessarily an integer. cp: Numeric vector of length tcp.max; the most probable locations of trend changepoints. These are derived by smoothing cpOccPr with a window of length tseg.min and picking up to prior$MaxKnotNum (subject to tcp.max) of the highest peaks in the smoothed curve. Entries may be NaN if insufficient changepoints are identified. Many entries may be false positives; they should not be treated as confirmed changepoints without further scrutiny. cpPr: Numeric vector of length tcp.max; the probabilities associated with the detected changepoints in cp. Remaining entries are NaN if ncp < tcp.max. cpCI: Numeric matrix of dimension tcp.max x 2; credible intervals for the detected changepoints in cp. cpAbruptChange: Numeric vector of length tcp.max; the jumps in the fitted trend at the detected changepoints. Y: Numeric vector of length N; posterior mean estimate of the trend component (Bayesian model average over sampled models). SD: Numeric vector of length N; posterior standard deviation of the trend estimate. CI: Numeric matrix of dimension N x 2; credible intervals (lower and upper envelopes) for the trend component. slp: Numeric vector of length N; time-varying slope of the fitted trend. slpSD: Numeric vector of length N; posterior standard deviation of the trend slope. slpSgnPosPr: Numeric vector of length N; probability that the trend slope is positive (i.e., increasing). For example, slpSgnPosPr[t] = 0.8 means that at time index t, 80% of sampled trend curves have positive slope. slpSgnZeroPr: Numeric vector of length N; probability that the trend slope is approximately zero (flat). The probability of a negative slope is given by 1 - slpSgnZeroPr - slpSgnPosPr. pos_ncp, neg_ncp, pos_ncpPr, neg_ncpPr, pos_cpOccPr, neg_cpOccPr, pos_cp, neg_cp, pos_cpPr, neg_cpPr, pos_cpAbruptChange, neg_cpAbruptChange, pos_cpCI, neg_cpCI: As above, but distinguishing changepoints where the trend jumps up (positive jump) versus jumps down (negative jump). For example, pos_ncp is the average number of changepoints with positive jumps. inc_ncp, dec_ncp, inc_ncpPr, dec_ncpPr, inc_cpOccPr, dec_cpOccPr, inc_cp, dec_cp, inc_cpPr, dec_cpPr, inc_cpAbruptChange, dec_cpAbruptChange, inc_cpCI, dec_cpCI: Analogous to the above, but distinguishing changepoints where the trend slope increases vs. decreases. For instance, if the slope changes from 0.4 to 2.5 at a changepoint, that changepoint contributes to inc_ncp.
season	List containing summaries for the estimated seasonal/periodic component (if present): ncp: Mean number of seasonal changepoints. ncpPr: Vector of length scp.max + 1, where scp.max = scp.minmax[2]; the probability of having 0, 1, ..., scp.max seasonal changepoints. cpOccPr: Numeric vector of length N; marginal probability that each time index is a seasonal changepoint. The same interpretation caveats as for the trend component apply. order: Numeric vector of length N; average harmonic order required to approximate the seasonal component. This is an average over sampled piecewise harmonic curves and is not necessarily an integer. cp: Numeric vector of length scp.max; the most probable seasonal changepoint locations, derived from the smoothed cpOccPr curve using a window size of sseg.min. If ncp < scp.max, the remaining entries are filled with NaN. cpPr: Numeric vector of length scp.max; probabilities associated with cp. Remaining entries are NaN if ncp < scp.max. cpCI: Numeric matrix of dimension scp.max x 2; credible intervals for the detected seasonal changepoints. cpAbruptChange: Numeric vector of length scp.max; jumps in the fitted seasonal component at the detected changepoints. Y: Numeric vector of length N; posterior mean estimate of the seasonal component. SD: Numeric vector of length N; posterior standard deviation of the seasonal estimate. CI: Numeric matrix of dimension N x 2; credible intervals (lower and upper envelopes) for the seasonal component. amp: Numeric vector of length N; time-varying amplitude of the seasonal component. ampSD: Numeric vector of length N; posterior standard deviation of the amplitude. pos_ncp, neg_ncp, pos_ncpPr, neg_ncpPr, pos_cpOccPr, neg_cpOccPr, pos_cp, neg_cp, pos_cpPr, neg_cpPr, pos_cpAbruptChange, neg_cpAbruptChange, pos_cpCI, neg_cpCI: Analogous to the trend component, but for seasonal changepoints that correspond to positive vs. negative jumps in the seasonal signal.
outlier	outlier is a list analogous to trend or season, but for the outlier component ( present only if setting hasOutlier=TRUE)

Note

The three functions beast(), beast.irreg(), and beast123() implement the same BEAST algorithm but expose different APIs:

• beast() operates on regular (evenly spaced) time series.

• beast.irreg() accepts irregular/unordered data and internally aggregates it to regular intervals.

• beast123() exposes a low-level interface via four lists: metadata, prior, mcmc, and extra.

There is a one-to-one correspondence between arguments of beast() / beast.irreg() and fields in metadata, prior, mcmc, and extra used by beast123(), Examples include:

start	<->	metadata$startTime
deltat	<->	metadata$deltaTime
deseasonalize	<->	metadata$deseasonalize
hasOutlier	<->	metadata$hasOutlier
scp.minmax[1:2]	<->	prior$seasonMinKnotNum, prior$seasonMaxKnotNum
sorder.minmax[1:2]	<->	prior$seasonMinOrder, prior$seasonMaxOrder
sseg.min	<->	prior$seasonMinSepDist
tcp.minmax[1:2]	<->	prior$trendMinKnotNum, prior$trendMaxKnotNum
torder.minmax[1:2]	<->	prior$trendMinOrder, prior$trendMaxOrder
tseg.leftmargin	<->	prior$trendLeftMargin
tseg.rightmargin	<->	prior$trendRightMargin
s.complexfct	<->	prior$seasonComplexityFactor
t.complexfct	<->	prior$trendComplexityFactor
mcmc.seed	<->	mcmc$seed
dump.ci	<->	extra$computeCredible.

Advanced users who need full control over all internal options are encouraged to use beast123() directly.

References

1. 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 change-point, trend, and seasonality in satellite time series data to track abrupt changes and nonlinear dynamics: A Bayesian ensemble algorithm. Remote Sensing of Environment, 232, 111181. (The main BEAST algorithm paper.)

2. 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, 102–119. (The Bayesian MCMC scheme used in BEAST.)

3. Hu, T., Toman, E. M., Chen, G., Shao, G., Zhou, Y., Li, Y., Zhao, K. and Feng, Y. (2021). Mapping fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, 250–261. (Example application of BEAST.)

See Also

beast, beast123, minesweeper, tetris, geeLandsat

Examples

library(Rbeast)



################################################################################
# Note: The BEAST algorithm is implemented for regular time series. 
# \code{beast.irreg} accepts irregular data but internally aggregates it to a 
# regular series before applying BEAST. For aggregation, both "time" and 
# "deltat" are needed to specify the original timestamps via "time" and the desired 
# regular interval via "deltat". If a cyclic component exists, "period" should also
# be provided; otherwise BEAST attempts to guess it via autocorrelation.
################################################################################

# 'ohio' is a data.frame containing an irregular Landsat time series of 
# reflectances and NDVI at an Ohio site. It includes multiple alternative 
# date formats: year (Y), month (M), day (D), day-of-year (doy), R "Date" 
# (rdate), and fractional year (time).

data(ohio)
str(ohio)
plot(ohio$rdate, ohio$ndvi, type = "o")  # NDVI is irregular and unordered in time

################################################################################
# Example 1: 'time' as numeric values
# 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, BEAST itself doesn't care about the unit. 
# So, the unit of 1/12 is irrelevant for BEAST. '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)

################################################################################
# Example 2: Aggregate to a monthly interval (deltat = 1/12) and provide 'period'
################################################################################

o <- beast.irreg(ohio$ndvi, time = ohio$time, deltat = 1/12, period = 1.0)
# Alternative (deprecated argument 'freq'):
# o <- beast.irreg(ohio$ndvi, time = ohio$time, deltat = 1/12, freq = 12)



## Not run: 

################################################################################
# Example 3: Aggregate at a half-monthly interval.
# Here period = 1: freq = period/deltat = 1/(1/24)=24 data points per period
################################################################################

o <- beast.irreg(ohio$ndvi, time = ohio$time, deltat = 1/24, period = 1)

################################################################################
# Example 4: 'time' as  R Dates (i.e, ohio$rdate). Unit is YEAR; 1/12 is ~1 month.
################################################################################

o <- beast.irreg(ohio$ndvi, time = ohio$rdate, deltat = 1/12)

################################################################################
# Example 5: 'time' as date strings. The unit is YEAR such that 1/12 is a month
################################################################################

o=beast.irreg(ohio$ndvi, time=ohio$datestr1,deltat=1/12) #"LT4-1984-03-27"  (YYYY-MM-DD)
o=beast.irreg(ohio$ndvi, time=ohio$datestr2,deltat=1/12) #"LT4-1984087ndvi" (YYYYDOY)
o=beast.irreg(ohio$ndvi, time=ohio$datestr3,deltat=1/12) #"1984,, 3/ 27"    (YYYY M D)

################################################################################
# Example 6: 'time' as date strings with explicit format specifiers
################################################################################

TIME <- list()
TIME$datestr <- ohio$datestr1
TIME$strfmt  <- "LT4-YYYY-MM-DD"    # e.g., "LT4-1984-03-27"
o <- beast.irreg(ohio$ndvi, time = TIME, deltat = 1/12)

TIME <- list()
TIME$datestr <- ohio$datestr2
TIME$strfmt  <- "LT4-YYYYDOYndvi"   # e.g., "LT4-1984087ndvi"
o <- beast.irreg(ohio$ndvi, time = TIME, deltat = 1/12)

################################################################################
# Example 7: 'time' as a list of date components
################################################################################

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)

Rbeast documentation built on Dec. 19, 2025, 9:06 a.m.