Bayesian changepoint detection detection and time series decomposition for trend, periodicity or seasonality, and abrupt changes
Description
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.
Usage
beast( y, start = 1, deltat = 1, 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'), detrend = FALSE, deseasonalize = FALSE, mcmc.seed = 0, mcmc.burnin = 200, mcmc.chains = 3, mcmc.thin = 5, mcmc.samples = 8000, ci = FALSE, precValue = 1.5, precPriorType = c('componentwise','uniform','constant','orderwise'), print.options = TRUE, print.progress = TRUE, quiet = FALSE, gui = FALSE, ... )beast( y, start = 1, deltat = 1, 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'), detrend = FALSE, deseasonalize = FALSE, mcmc.seed = 0, mcmc.burnin = 200, mcmc.chains = 3, mcmc.thin = 5, mcmc.samples = 8000, ci = FALSE, precValue = 1.5, precPriorType = c('componentwise','uniform','constant','orderwise'), print.options = TRUE, print.progress = TRUE, quiet = FALSE, gui = FALSE, ... )
Arguments
y |
a vector for an evenly-spaced regular time series. Missing values such as NA and NaN are allowed.
If a list of multiple time series is provided for |
start |
numeric (default to 1.0) or |
deltat |
numeric (default to 1.0) or string; the time interval between consecutive data points. Its unit should be consistent with |
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 0-th 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 burn-in 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 |
ci |
boolean; If |
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.
|
print.options |
boolean. If |
print.progress |
boolean;If |
quiet |
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 |
Value
The output is an object of class "beast". It is a list, consisting of the following variables. Its structure is the same as the outputs from the other two alternative functions beast.irreg and beast123. 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 greater marg_lik, the better the fitting for a given time series; that is, -10 will be better than -100 and 10 better than any of them. |
R2 |
numeric; the R-square 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:
|
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, minesweeper, tetris, geeLandsat
Examples
library(Rbeast) #------------------------------------Example 1----------------------------------------# # 'googletrend_beach' is the monthly Google Trend popularity of searching for 'beach' # in the US from 2004 to 2022. Sudden changes in the time series coincide with known # extreme weather events (e.g., 2006 North American Blizzard, 2011 US hottest summer # on record, Record warm January in 2016) or the covid19 outbreak. out <- beast(googletrend_beach) plot(out) plot(out, vars=c('t','slpsgn') ) # plot the trend and probability of slope sign only. # In the slpsgn panel, the upper red portion refers to # probability of trend slope being positive, the middle # green to the prob of slope being zero, and the lower # blue to the probability of slope being negative. # Run "?plot.beast" for details on the plot function. #------------------------------------Example 2----------------------------------------# # Yellowstone is a half-monthly satellite time series of 774 NDVI(vegetation greeness) # observations starting from July 1-15,1981(i.e., start=c(1981,7,7)) at a Yellowstone # forest site. It has 24 data points per year (i.e., freq=24). Note that the beast # function hanldes only evenly-spaced regular time series. Irregular data need to be # first aggegrated at a regular time interval of your choice--the aggregation # functionality is implemented in beast.irreg() and beast123(). data(Yellowstone) plot( 1981.5+(0:773)/24, Yellowstone, type='l') # A sudden drop in greenness in 1988 # due to the 1988 Yellowstone Fire # Yellowstone is not a object of class 'ts' but a pure vector without time attributes. # Below, no extra argument is supplied, so default values (i.e.,start=1, deltat=1) are # used and the time is 1:774. 'period' is missing and so is guessed via auto-correlation. # Use of auto-correlation to compute the period of a cyclic time series is not always # reliable, so it is suggested to always supply 'period' directly, as in Example 2 and # Example 3. o = beast(Yellowstone) # By defualt, the times assumed to be 1:length(Yellowstone) # and a periodic component is assumed (season='harmonic') plot(o) #o = beast(Yellowstone, quiet=TRUE) # print no warning messages #o = beast(Yellowstone, quiet=TRUE, print.progress=FALSE) # print nothing #------------------------------------Example 3----------------------------------------# # The time info such as start,delta,and period is explicitly provided. 'start' can be # given as (1) a fractional number, (2) a vector comprising year, month,& day, or (3) # a R's Date. In (1), the unit of start and deltat does not necessarily refer to time and can # be arbitrary (e.g., a sequence of data observed at evenly-spaced distaces along a # transect or a elevation gradient) # (1) Unknown unit such that 1981.5137 can be interpreted arbitrarily o=beast(Yellowstone, start=1981.5137, deltat=1/24, period=1.0) # Use a string to explictly specify a time unit so that times are intepreted as dates # o=beast(Yellowstone, start=1981.5137, deltat='1/24 year', period=1.0) # 1.0 = 1 yr # o=beast(Yellowstone, start=1981.5137, deltat='0.5 mon', period=1.0) # 1.0 = 1 yr # o=beast(Yellowstone, start=1981.5137, deltat=1/24, period='1 yr') # 1/24 = 1/24 yr # o=beast(Yellowstone, start=1981.5137, deltat=1/24, period='365 days')# 1/24 = 1/24 yr # (2) start is provided as YMD, the unit is year: deltat=1/24 year=0.5 month # o=beast(Yellowstone, start=c(1981,7,7), deltat=1/24, period=1.0) # (3) start is provided as Date, the unit is year: deltat=1/24 year=0.5 month #o=beast(Yellowstone, start=as.Date('1981-7-7'), deltat=1/24, period=1.0) print(o) # o is a R LIST object with many fields str(o) # See a list of fields in o plot(o) # plot many variables plot(o, vars=c('y','s','t') ) # plot the Y, seasonal, and trend components only plot(o, vars=c('s','scp','samp','t','tcp','tslp'))# Plot some selected variables in # 'o'. Type "?plot.beast" to see # more about vars plot(o, vars=c('s','t'),col=c('red','blue') ) # Specify colors of selected subplots plot(o$time, o$season$Y,type='l') # directly plot output: the fitted season plot(o$time, o$season$cpOccPr) # directly plot output: season chgpt prob plot(o$time, o$trend$Y,type='l') # directly plot output: the fitted trend plot(o$time, o$trend$cpOccPr) # directly plot output: trend chgpt occurrence prob plot(o$time, o$season$order) # directly plot output: avg harmonic order plot(o, interactive=TRUE) # manually choose which variables to plot #------------------------------------Example 4----------------------------------------# # Specify other arguments explicitly. Default values are used for missing parameters. # Note that beast(), beast.irreg(), and beast123() call the same internal C/C++ library, # so in beast(), the input parameters will be converted to metadata, prior, mcmc, and # extra parameters as explained for the beast123() function. Or type 'View(beast)' to # check the parameter assignment in the code. # In R's terminology, the number of datapoints per period is also called 'freq'. In this # version, the 'freq' argument is obsolete and replaced by 'period'. # period=deltat*number_of_datapoints_per_period = 1.0*24=24 because deltat is set to 1.0 by # default and this signal has 24 samples per period. out = beast(Yellowstone, period=24.0, mcmc.samples=5000, tseg.min=20) # period=deltat*number_of_datapoints_per_period = 1/24*24=1.0. # out = beast(Yellowstone, deltat=1/24 period=1.0, mcmc.samples=5000, tseg.min=20) out = beast( Yellowstone, # Yellowstone: a pure numeric vector wo time info start = 1981.51, deltat = 1/24, period = 1.0, # period=delta*number_of_datapoints_per_period season = 'harmonic', # periodic compnt exisits,fitted as a harmonic curve scp.minmax = c(0,3), # min and max numbers of seasonal changpts allowed sorder.minmax = c(1,5), # min and max harmonic orders allowed sseg.min = 24, # the min length of segments btw neighboring chnpts # '24' means 24 datapoints; the unit is datapoint. sseg.leftmargin= 40, # no seasonal chgpts allowed in the starting 40 datapoints tcp.minmax = c(0,10),# min and max numbers of changpts allowed in the trend torder.minmax = c(0,1), # min and maxx polynomial orders to fit trend tseg.min = 24, # the min length of segments btw neighboring trend chnpts tseg.leftmargin= 10, # no trend chgpts allowed in the starting 10 datapoints deseasonalize = TRUE, # remove the global seasonality before fitting the beast model detrend = TRUE, # remove the global trend before fitting the beast model mcmc.seed = 0, # a seed for mcmc's random nummber generator; use a # non-zero integer to reproduce results across runs mcmc.burnin = 500, # number of initial iterations discarded mcmc.chains = 2, # number of chains mcmc.thin = 3, # include samples every 3 iterations mcmc.samples = 6000 # number of samples taken per chain # total iteration: (500+3*6000)*2 ) plot(out) plot(out, interactive=TRUE) #------------------------------------Example 5----------------------------------------# # Run an interactive GUI to visualize how BEAST is samplinig from the possible model # spaces in terms of the numbers and timings of seasonal and trend changepoints. # The GUI inferface allows changing the option parameters interactively. This GUI is # only available on Win x64 machines, not Mac or Linux. ## Not run: beast(Yellowstone, period=24, gui=TRUE) ## End(Not run) #------------------------------------Example 6----------------------------------------# # Apply beast to trend-only data. 'Nile' is the ANNUAL river flow of the river # Nile at Aswan since 1871. It is a 'ts' object; its time attributes (start=1871, # end=1970,frequency=1) are used to replace the user-supplied start,deltat, and freq, # if any. data(Nile) plot(Nile) attributes(Nile) # a ts object with time attributes (i.e., tsp=(start,end,freq) o = beast(Nile) # start=1871, delta=1, and freq=1 taken from Nile itself plot(o) o = beast(Nile, # the same as above. The user-supplied values (i.e., 2023, start=2023, # 9999) are ignored bcz Nile carries its own time attributes. period=9999, # Its frequency tag is 1 (i.e., trend-only), so season='none' season='harmonic' # is used instead of the supplied 'harmonic' ) #------------------------------------Example 7----------------------------------------# # NileVec is a pure data vector. The first run below is WRONG bcz NileVec was assumed # to have a perodic component by default and beast gets a best estimate of freq=6 while # the true value is freq=1. To fit a trend-only model, season='none' has to be explicitly # specified, as in the 2nd & 3rd funs. NileVec = as.vector(Nile) # NileVec is not a ts obj but a pure numeric data vector o = beast(NileVec) # WRONG WAY to call: No time attributes available to interpret # NileVec. By default, beast assumes season='harmonic', start=1, # & deltat=1. 'freq' is missing and guessed to be 6 (WRONG). plot(o) # WRONG Results: The result has a suprious seasonal component o=beast(NileVec,season='none') # The correct way to call: Use season='none' for trend-only # analysis; the default time is the integer indices # "1:length(NileVec)'. print(o$time) o=beast(NileVec, # Recommended way to call: The true time attributes are start = 1871, # given explicitly through start and deltat (or freq if deltat = 1, # there is a cyclic/seasonal cmponent). season = 'none') print(o$time) plot(o) #------------------------------------Example 8----------------------------------------# # beast can handle missing data. co2 is a monthly time series (i.e.,freq=12) starting # from Jan 1959. We generate some missing values at random indices ## Not run: data(co2) attributes(co2) # A ts object with time attributes (i.e., tsp) badIdx = sample( 1:length(co2), 50) # Get a set of random indices co2[badIdx] = NA # Insert some data gaps out=beast(co2) # co2 is a ts object and its 'tsp' time attributes are used to get the # true time info. No need to specify 'start','deltat', & freq explicity. out=beast(co2, # The supplied time/period values will be ignored bcz start = c(1959,1,15),# co2 is a ts object; the correct period = 1 will be deltat = 1/12, # used. period = 365) print(out) plot(out) ## End(Not run) #------------------------------------Example 9----------------------------------------# # Apply beast to time seris-like sequence data: the unit of sequences is not # necessarily time. data(CNAchrom11) # DNA copy number alterations in Chromesome 11 for cell line GM05296 # The data is orderd by genomic position (not time), and the values # are the log2-based intensity ratio of copy numbers between the sample # the reference. A value of zero means no gain or loss in copy number. o = beast(CNAchrom11,season='none') # season is a misnomer here bcz the data has nothing # to do with time. Regardless, we fit only a trend. plot(o) #------------------------------------Example 10---------------------------------------# # Apply beast to time seris-like data: the unit of sequences is not necessarily time. # Age of Death of Successive Kings of England # If the data link is deprecated, install the time series data library instead, # which is available at https://pkg.yangzhuoranyang.com/tsdl/ # install.packages("devtools") # devtools::install_github("FinYang/tsdl") # kings = tsdl::tsdl[[293]] kings = scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3) out = beast(kings,season='none') plot(out) #------------------------------------Example 11---------------------------------------# # Another example from the tsdl data library # Number of monthly births in New York from Jan 1946 to Dec 1959 # If the data link becomes invalid, install the time series data package instead # install.packages("devtools") # devtools::install_github("FinYang/tsdl") # kings = tsdl::tsdl[[534]] births = scan("http://robjhyndman.com/tsdldata/data/nybirths.dat") out = beast(births,start=c(1946,1,15), deltat=1/12 ) plot(out) # the result is wrong bcz the guessed freq via auto-correlation by beast # is 2 rather than 12, so we recommend always specifying 'freq' explicitly # for those time series with a periodic component, as shown below. out = beast(births,start=c(1946,1,15), deltat=1/12, freq =12 ) out = beast(births,start=c(1946,1,15), deltat=1/12, period=1.0 ) plot(out) #------------------------------------Example 12---------------------------------------# # Daily confirmed COVID-19 new cases and deaths across the globe ## Not run: data(covid19) plot(covid19$date, covid19$newcases, type='l') newcases = sqrt( covid19$newcases ) # Apply a square root-transformation # This ts varies periodically every 7 days. 7 days can't be precisely represented # in the unit of year bcz some years has 365 days and others has 366. BEAST can hanlde # this in two ways. #(1) Use the date number as the time unit--the num of days lapsed since 1970-01-01. datenum = as.numeric(covid19$date) o = beast(newcases, start=min(datenum), deltat=1, period=7) o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date. plot(o) #(2) Use strings to explicitly specify deltat and period with a unit. startdate = covid19$date[1] o = beast(newcases, start=startdate, deltat='1day', period='7days') plot(o) ## End(Not run) #------------------------------------Example 13---------------------------------------# # The old API interface of beast is still made available but NOT recommended. It is # kept mainly to ensure the working of the sample code on Page 475 in the text # Ecological Metods by Drs. Southwood and Henderson. ## Not run: # The interface as shown here will be deprecated and NOT recommended. beast(Yellowstone, 24) #24 is the freq: number of datapoints per period # Specify the model or MCMC parameters through opt as in Rbeast v0.2 opt=list() #Create an empty list to append individual model parameters opt$period=24 #Period of the cyclic component (i.e.,freq in the new version) opt$minSeasonOrder=2 #Min harmonic order allowed in fitting season component opt$maxSeasonOrder=8 #Max harmonic order allowed in fititing season component opt$minTrendOrder=0 #Min polynomial order allowed to fit trend (0 for constant) opt$maxTrendOrder=1 #Max polynomial order allowed to fit trend (1 for linear term) opt$minSepDist_Season=20#Min separation time btw neighboring season changepoints opt$minSepDist_Trend=20 #Min separation time btw neighboring trend changepoints opt$maxKnotNum_Season=4 #Max number of season changepoints allowed opt$maxKnotNum_Trend=10 #Max number of trend changepoints allowed opt$omittedValue=NA #A customized value to indicate bad/missing values in the time #series, in additon to those NA or NaN values. # The following parameters used to configure the reverisible-jump MCMC (RJMCC) sampler opt$chainNumber=2 #Number of parallel MCMC chains opt$sample=1000 #Number of samples to be collected per chain opt$thinningFactor=3 #A factor to thin chains opt$burnin=500 #Number of burn-in samples discarded at the start opt$maxMoveStepSize=30 #For the move proposal, the max window allowed in jumping from #the current changepoint opt$resamplingSeasonOrderProb=0.2 #The probability of selecting a re-sampling proposal #(e.g., resample seasonal harmonic order) opt$resamplingTrendOrderProb=0.2 #The probability of selecting a re-sampling proposal #(e.g., resample trend polynomial order) opt$seed=65654 #A seed for the random generator: If seed=0,random numbers differ #for different BEAST runs. Setting seed to a chosen non-zero integer #will allow reproducing the same result for different BEAST runs. beast(Yellowstone, opt) ## End(Not run)library(Rbeast) #------------------------------------Example 1----------------------------------------# # 'googletrend_beach' is the monthly Google Trend popularity of searching for 'beach' # in the US from 2004 to 2022. Sudden changes in the time series coincide with known # extreme weather events (e.g., 2006 North American Blizzard, 2011 US hottest summer # on record, Record warm January in 2016) or the covid19 outbreak. out <- beast(googletrend_beach) plot(out) plot(out, vars=c('t','slpsgn') ) # plot the trend and probability of slope sign only. # In the slpsgn panel, the upper red portion refers to # probability of trend slope being positive, the middle # green to the prob of slope being zero, and the lower # blue to the probability of slope being negative. # Run "?plot.beast" for details on the plot function. #------------------------------------Example 2----------------------------------------# # Yellowstone is a half-monthly satellite time series of 774 NDVI(vegetation greeness) # observations starting from July 1-15,1981(i.e., start=c(1981,7,7)) at a Yellowstone # forest site. It has 24 data points per year (i.e., freq=24). Note that the beast # function hanldes only evenly-spaced regular time series. Irregular data need to be # first aggegrated at a regular time interval of your choice--the aggregation # functionality is implemented in beast.irreg() and beast123(). data(Yellowstone) plot( 1981.5+(0:773)/24, Yellowstone, type='l') # A sudden drop in greenness in 1988 # due to the 1988 Yellowstone Fire # Yellowstone is not a object of class 'ts' but a pure vector without time attributes. # Below, no extra argument is supplied, so default values (i.e.,start=1, deltat=1) are # used and the time is 1:774. 'period' is missing and so is guessed via auto-correlation. # Use of auto-correlation to compute the period of a cyclic time series is not always # reliable, so it is suggested to always supply 'period' directly, as in Example 2 and # Example 3. o = beast(Yellowstone) # By defualt, the times assumed to be 1:length(Yellowstone) # and a periodic component is assumed (season='harmonic') plot(o) #o = beast(Yellowstone, quiet=TRUE) # print no warning messages #o = beast(Yellowstone, quiet=TRUE, print.progress=FALSE) # print nothing #------------------------------------Example 3----------------------------------------# # The time info such as start,delta,and period is explicitly provided. 'start' can be # given as (1) a fractional number, (2) a vector comprising year, month,& day, or (3) # a R's Date. In (1), the unit of start and deltat does not necessarily refer to time and can # be arbitrary (e.g., a sequence of data observed at evenly-spaced distaces along a # transect or a elevation gradient) # (1) Unknown unit such that 1981.5137 can be interpreted arbitrarily o=beast(Yellowstone, start=1981.5137, deltat=1/24, period=1.0) # Use a string to explictly specify a time unit so that times are intepreted as dates # o=beast(Yellowstone, start=1981.5137, deltat='1/24 year', period=1.0) # 1.0 = 1 yr # o=beast(Yellowstone, start=1981.5137, deltat='0.5 mon', period=1.0) # 1.0 = 1 yr # o=beast(Yellowstone, start=1981.5137, deltat=1/24, period='1 yr') # 1/24 = 1/24 yr # o=beast(Yellowstone, start=1981.5137, deltat=1/24, period='365 days')# 1/24 = 1/24 yr # (2) start is provided as YMD, the unit is year: deltat=1/24 year=0.5 month # o=beast(Yellowstone, start=c(1981,7,7), deltat=1/24, period=1.0) # (3) start is provided as Date, the unit is year: deltat=1/24 year=0.5 month #o=beast(Yellowstone, start=as.Date('1981-7-7'), deltat=1/24, period=1.0) print(o) # o is a R LIST object with many fields str(o) # See a list of fields in o plot(o) # plot many variables plot(o, vars=c('y','s','t') ) # plot the Y, seasonal, and trend components only plot(o, vars=c('s','scp','samp','t','tcp','tslp'))# Plot some selected variables in # 'o'. Type "?plot.beast" to see # more about vars plot(o, vars=c('s','t'),col=c('red','blue') ) # Specify colors of selected subplots plot(o$time, o$season$Y,type='l') # directly plot output: the fitted season plot(o$time, o$season$cpOccPr) # directly plot output: season chgpt prob plot(o$time, o$trend$Y,type='l') # directly plot output: the fitted trend plot(o$time, o$trend$cpOccPr) # directly plot output: trend chgpt occurrence prob plot(o$time, o$season$order) # directly plot output: avg harmonic order plot(o, interactive=TRUE) # manually choose which variables to plot #------------------------------------Example 4----------------------------------------# # Specify other arguments explicitly. Default values are used for missing parameters. # Note that beast(), beast.irreg(), and beast123() call the same internal C/C++ library, # so in beast(), the input parameters will be converted to metadata, prior, mcmc, and # extra parameters as explained for the beast123() function. Or type 'View(beast)' to # check the parameter assignment in the code. # In R's terminology, the number of datapoints per period is also called 'freq'. In this # version, the 'freq' argument is obsolete and replaced by 'period'. # period=deltat*number_of_datapoints_per_period = 1.0*24=24 because deltat is set to 1.0 by # default and this signal has 24 samples per period. out = beast(Yellowstone, period=24.0, mcmc.samples=5000, tseg.min=20) # period=deltat*number_of_datapoints_per_period = 1/24*24=1.0. # out = beast(Yellowstone, deltat=1/24 period=1.0, mcmc.samples=5000, tseg.min=20) out = beast( Yellowstone, # Yellowstone: a pure numeric vector wo time info start = 1981.51, deltat = 1/24, period = 1.0, # period=delta*number_of_datapoints_per_period season = 'harmonic', # periodic compnt exisits,fitted as a harmonic curve scp.minmax = c(0,3), # min and max numbers of seasonal changpts allowed sorder.minmax = c(1,5), # min and max harmonic orders allowed sseg.min = 24, # the min length of segments btw neighboring chnpts # '24' means 24 datapoints; the unit is datapoint. sseg.leftmargin= 40, # no seasonal chgpts allowed in the starting 40 datapoints tcp.minmax = c(0,10),# min and max numbers of changpts allowed in the trend torder.minmax = c(0,1), # min and maxx polynomial orders to fit trend tseg.min = 24, # the min length of segments btw neighboring trend chnpts tseg.leftmargin= 10, # no trend chgpts allowed in the starting 10 datapoints deseasonalize = TRUE, # remove the global seasonality before fitting the beast model detrend = TRUE, # remove the global trend before fitting the beast model mcmc.seed = 0, # a seed for mcmc's random nummber generator; use a # non-zero integer to reproduce results across runs mcmc.burnin = 500, # number of initial iterations discarded mcmc.chains = 2, # number of chains mcmc.thin = 3, # include samples every 3 iterations mcmc.samples = 6000 # number of samples taken per chain # total iteration: (500+3*6000)*2 ) plot(out) plot(out, interactive=TRUE) #------------------------------------Example 5----------------------------------------# # Run an interactive GUI to visualize how BEAST is samplinig from the possible model # spaces in terms of the numbers and timings of seasonal and trend changepoints. # The GUI inferface allows changing the option parameters interactively. This GUI is # only available on Win x64 machines, not Mac or Linux. ## Not run: beast(Yellowstone, period=24, gui=TRUE) ## End(Not run) #------------------------------------Example 6----------------------------------------# # Apply beast to trend-only data. 'Nile' is the ANNUAL river flow of the river # Nile at Aswan since 1871. It is a 'ts' object; its time attributes (start=1871, # end=1970,frequency=1) are used to replace the user-supplied start,deltat, and freq, # if any. data(Nile) plot(Nile) attributes(Nile) # a ts object with time attributes (i.e., tsp=(start,end,freq) o = beast(Nile) # start=1871, delta=1, and freq=1 taken from Nile itself plot(o) o = beast(Nile, # the same as above. The user-supplied values (i.e., 2023, start=2023, # 9999) are ignored bcz Nile carries its own time attributes. period=9999, # Its frequency tag is 1 (i.e., trend-only), so season='none' season='harmonic' # is used instead of the supplied 'harmonic' ) #------------------------------------Example 7----------------------------------------# # NileVec is a pure data vector. The first run below is WRONG bcz NileVec was assumed # to have a perodic component by default and beast gets a best estimate of freq=6 while # the true value is freq=1. To fit a trend-only model, season='none' has to be explicitly # specified, as in the 2nd & 3rd funs. NileVec = as.vector(Nile) # NileVec is not a ts obj but a pure numeric data vector o = beast(NileVec) # WRONG WAY to call: No time attributes available to interpret # NileVec. By default, beast assumes season='harmonic', start=1, # & deltat=1. 'freq' is missing and guessed to be 6 (WRONG). plot(o) # WRONG Results: The result has a suprious seasonal component o=beast(NileVec,season='none') # The correct way to call: Use season='none' for trend-only # analysis; the default time is the integer indices # "1:length(NileVec)'. print(o$time) o=beast(NileVec, # Recommended way to call: The true time attributes are start = 1871, # given explicitly through start and deltat (or freq if deltat = 1, # there is a cyclic/seasonal cmponent). season = 'none') print(o$time) plot(o) #------------------------------------Example 8----------------------------------------# # beast can handle missing data. co2 is a monthly time series (i.e.,freq=12) starting # from Jan 1959. We generate some missing values at random indices ## Not run: data(co2) attributes(co2) # A ts object with time attributes (i.e., tsp) badIdx = sample( 1:length(co2), 50) # Get a set of random indices co2[badIdx] = NA # Insert some data gaps out=beast(co2) # co2 is a ts object and its 'tsp' time attributes are used to get the # true time info. No need to specify 'start','deltat', & freq explicity. out=beast(co2, # The supplied time/period values will be ignored bcz start = c(1959,1,15),# co2 is a ts object; the correct period = 1 will be deltat = 1/12, # used. period = 365) print(out) plot(out) ## End(Not run) #------------------------------------Example 9----------------------------------------# # Apply beast to time seris-like sequence data: the unit of sequences is not # necessarily time. data(CNAchrom11) # DNA copy number alterations in Chromesome 11 for cell line GM05296 # The data is orderd by genomic position (not time), and the values # are the log2-based intensity ratio of copy numbers between the sample # the reference. A value of zero means no gain or loss in copy number. o = beast(CNAchrom11,season='none') # season is a misnomer here bcz the data has nothing # to do with time. Regardless, we fit only a trend. plot(o) #------------------------------------Example 10---------------------------------------# # Apply beast to time seris-like data: the unit of sequences is not necessarily time. # Age of Death of Successive Kings of England # If the data link is deprecated, install the time series data library instead, # which is available at https://pkg.yangzhuoranyang.com/tsdl/ # install.packages("devtools") # devtools::install_github("FinYang/tsdl") # kings = tsdl::tsdl[[293]] kings = scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3) out = beast(kings,season='none') plot(out) #------------------------------------Example 11---------------------------------------# # Another example from the tsdl data library # Number of monthly births in New York from Jan 1946 to Dec 1959 # If the data link becomes invalid, install the time series data package instead # install.packages("devtools") # devtools::install_github("FinYang/tsdl") # kings = tsdl::tsdl[[534]] births = scan("http://robjhyndman.com/tsdldata/data/nybirths.dat") out = beast(births,start=c(1946,1,15), deltat=1/12 ) plot(out) # the result is wrong bcz the guessed freq via auto-correlation by beast # is 2 rather than 12, so we recommend always specifying 'freq' explicitly # for those time series with a periodic component, as shown below. out = beast(births,start=c(1946,1,15), deltat=1/12, freq =12 ) out = beast(births,start=c(1946,1,15), deltat=1/12, period=1.0 ) plot(out) #------------------------------------Example 12---------------------------------------# # Daily confirmed COVID-19 new cases and deaths across the globe ## Not run: data(covid19) plot(covid19$date, covid19$newcases, type='l') newcases = sqrt( covid19$newcases ) # Apply a square root-transformation # This ts varies periodically every 7 days. 7 days can't be precisely represented # in the unit of year bcz some years has 365 days and others has 366. BEAST can hanlde # this in two ways. #(1) Use the date number as the time unit--the num of days lapsed since 1970-01-01. datenum = as.numeric(covid19$date) o = beast(newcases, start=min(datenum), deltat=1, period=7) o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date. plot(o) #(2) Use strings to explicitly specify deltat and period with a unit. startdate = covid19$date[1] o = beast(newcases, start=startdate, deltat='1day', period='7days') plot(o) ## End(Not run) #------------------------------------Example 13---------------------------------------# # The old API interface of beast is still made available but NOT recommended. It is # kept mainly to ensure the working of the sample code on Page 475 in the text # Ecological Metods by Drs. Southwood and Henderson. ## Not run: # The interface as shown here will be deprecated and NOT recommended. beast(Yellowstone, 24) #24 is the freq: number of datapoints per period # Specify the model or MCMC parameters through opt as in Rbeast v0.2 opt=list() #Create an empty list to append individual model parameters opt$period=24 #Period of the cyclic component (i.e.,freq in the new version) opt$minSeasonOrder=2 #Min harmonic order allowed in fitting season component opt$maxSeasonOrder=8 #Max harmonic order allowed in fititing season component opt$minTrendOrder=0 #Min polynomial order allowed to fit trend (0 for constant) opt$maxTrendOrder=1 #Max polynomial order allowed to fit trend (1 for linear term) opt$minSepDist_Season=20#Min separation time btw neighboring season changepoints opt$minSepDist_Trend=20 #Min separation time btw neighboring trend changepoints opt$maxKnotNum_Season=4 #Max number of season changepoints allowed opt$maxKnotNum_Trend=10 #Max number of trend changepoints allowed opt$omittedValue=NA #A customized value to indicate bad/missing values in the time #series, in additon to those NA or NaN values. # The following parameters used to configure the reverisible-jump MCMC (RJMCC) sampler opt$chainNumber=2 #Number of parallel MCMC chains opt$sample=1000 #Number of samples to be collected per chain opt$thinningFactor=3 #A factor to thin chains opt$burnin=500 #Number of burn-in samples discarded at the start opt$maxMoveStepSize=30 #For the move proposal, the max window allowed in jumping from #the current changepoint opt$resamplingSeasonOrderProb=0.2 #The probability of selecting a re-sampling proposal #(e.g., resample seasonal harmonic order) opt$resamplingTrendOrderProb=0.2 #The probability of selecting a re-sampling proposal #(e.g., resample trend polynomial order) opt$seed=65654 #A seed for the random generator: If seed=0,random numbers differ #for different BEAST runs. Setting seed to a chosen non-zero integer #will allow reproducing the same result for different BEAST runs. beast(Yellowstone, opt) ## End(Not run)
Bayesian time series decomposition for changepoint, trend, and periodicity or seasonality
Description
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.
Usage
beast.irreg(\ y, time, deltat = NULL, period = NULL, season = c("harmonic", "svd", "dummy", "none"), 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'), detrend = FALSE, deseasonalize = FALSE, mcmc.seed = 0, mcmc.burnin = 200, mcmc.chains = 3, mcmc.thin = 5, mcmc.samples = 8000, ci = FALSE, precValue = 1.5, precPriorType = c('componentwise','uniform','constant','orderwise'), print.options = TRUE, print.progress = TRUE, quiet = FALSE, gui = FALSE, ... )beast.irreg(\ y, time, deltat = NULL, period = NULL, season = c("harmonic", "svd", "dummy", "none"), 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'), detrend = FALSE, deseasonalize = FALSE, mcmc.seed = 0, mcmc.burnin = 200, mcmc.chains = 3, mcmc.thin = 5, mcmc.samples = 8000, ci = FALSE, precValue = 1.5, precPriorType = c('componentwise','uniform','constant','orderwise'), print.options = TRUE, print.progress = TRUE, quiet = FALSE, gui = FALSE, ... )
Arguments
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 |
period |
numeric or string. Specify the period for the periodic/seasonal component in |
season |
characters (default to 'harmonic'); specify if
|
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. |
torder.minmax |
a vector of 2 integers (>=0); the min and max orders of the polynomials considered to fit the trend component. The 0-th order corresponds to a constant term/a flat line and the 1st order is a line. If |
tcp.minmax |
a vector of 2 integers; 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]. |
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 burn-in 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 |
ci |
boolean; If |
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.
|
print.options |
boolean. If |
print.progress |
boolean;If |
quiet |
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 |
Value
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 R-square 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:
|
Note
x
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast123, minesweeper, tetris, geeLandsat
Examples
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 auto-correlation ###################################################################################### # '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 half-monthly 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) #"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) ###################################################################################### # 'time' is given as data strings, with a format specifier TIME =list() TIME$datestr = ohio$datestr1 TIME$strfmt = "LT4-YYYY-MM-DD" # "LT4-1984-03-27" o=beast.irreg(ohio$ndvi, time=TIME,deltat=1/12) TIME =list() TIME$datestr = ohio$datestr2 TIME$strfmt = "LT4-YYYYDOYndvi" # LT4-1984087ndvi 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)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 auto-correlation ###################################################################################### # '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 half-monthly 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) #"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) ###################################################################################### # 'time' is given as data strings, with a format specifier TIME =list() TIME$datestr = ohio$datestr1 TIME$strfmt = "LT4-YYYY-MM-DD" # "LT4-1984-03-27" o=beast.irreg(ohio$ndvi, time=TIME,deltat=1/12) TIME =list() TIME$datestr = ohio$datestr2 TIME$strfmt = "LT4-YYYYDOYndvi" # LT4-1984087ndvi 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)
Bayesian time series decomposition for changepoint, trend, and periodicity or seasonality
Description
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.
Usage
beast123( Y, metadata=list(), prior =list(), mcmc =list(), extra =list(), season = c('harmonic','svd','dummy','none'), method = c( 'bayes', 'bic', 'aic', 'aicc','hic'), ...)beast123( Y, metadata=list(), prior =list(), mcmc =list(), extra =list(), season = c('harmonic','svd','dummy','none'), method = c( 'bayes', 'bic', 'aic', 'aicc','hic'), ...)
Arguments
Y |
a 1D vector, 2D matrix, or 3D array of numeric data. Missing values are allowed and can be indicated by
|
metadata |
(optional). If present,
|
prior |
(optional). a list object consisting of the hyperprior parameters in the Bayesian formulation of the BEAST model. Because they are part of the model, the fitting result may be sensitive to the choices of these hyperparameters. If
|
mcmc |
(optional). a list object consisting of parameters to configure the MCMC inference. These parameter are not part of the Bayesian formulation of the BEAST model but are the settings for the reversible-jump MCMC to generate MCMC chains. Due to the MCMC nature, the longer the simulation chain is, the better the fitting result. Below are possible parameters:
|
extra |
(optional). a list object consisting of flags to control the outputs from the BEAST runs or configure other program setting. Below are possible parameters:
|
season |
characters (default to 'harmonic'); specify if
|
method |
a string (default to 'bayes'); specify the method for formulating model posterior probability.
|
... |
additional parameters, not used currently but reserved for future extension |
Value
The output is an object of class "beast". It is a list, consisting of the following variables. Exact sizes of the variables depend on the types of the input Y as well as the specified output time dimension extra$whichOutputDimIsTime. In the explanations below, we assume the input Y is a single time series of length N; the dimensions for 2D or 2D inputs may be interpreted accordingly:
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 R-square 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 numeric consisting of various outputs related to the estimated trend component:
|
season |
a list object numeric consisting of various outputs related to the estimated seasonal/periodic component:
|
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, minesweeper, tetris, geeLandsat
Examples
#--------------------------------NOTE-------------------------------------------------# # beast123() is an all-inclusive function that duplicates the functionalities of beast # and beast.irreg. It can handle a single, multiple, or 3D of stacked time series, being # either regular or irregular. It allows for customization through four LIST arguments: # metadata -- additional info about the input Y # prior -- prior parameters for the beast model # mcmc -- MCMC simulation setting # extra -- misc parameters turning on/off outputs and setting up parallel computations # # Despite being essentially the same as beast and beast.irreg, beast123 is provided mainly # to support concurrent handling of multiple time series (e.g., stacked satellite images) # via parallel computing: When processing stacked raster layers, DO NOT iterate pixel by pixel # using beast() or beast.irreg() via an external parallel caller (e.g., doParallel or foreach). # Instread, please use beast123(), which supports mulithreading internally. #------------------------------Example 1: one time series with seasonalty-------------# # Yellowstone is a half-monthly time series of 774 NDVI measurmments at a Yellowstone # site starting from July 1-15,1981(i.e., start=c(1981,7,7). It has 24 data points per # year (freq=24). library(Rbeast) data(Yellowstone) plot(Yellowstone) # Below, the four option args are missing, so defalut values will be used, with some # warning messages given to altert this. By default, the input Y is assumed to be regular # with a seasonal component. The default arg values used will be printed out and they can # serve as a template to customize the parameters. o = beast123(Yellowstone) plot(o) #------------------------------Example 2: a trend-only time series-------------------# # Nile is an annual river flow time series (i.e., no periodic variation). So, season # is set to 'none' to indicate trend-only analysis. Default values are used for other # missing options. Unlike the beast() function, beast123 does NOT use the time attributes # of a 'ts' object. For example, Nile is treated as a pure data number; its (start=1871, # end=1970, freq=1) attributes are ignored. The default times 1:length(Nile) are used # instead. The true time info need to be specified by the 'metadata' parameter, as shown # in the next example. o = beast123(Nile,season='none') plot(o) #------------------------------Example 3: call via the full API interface-----------# # Specify metadata, prior, mcmc, and extra explicitly. Only 'prior' is the true statistical # model parameters of BEAST; the other three are just options to configure the input/ouput # or the computation process. ## Not run: # metadata is NOT part of BEAST itself, but some extra info to describe the input # time series Y. Below, the input Y is the 'Yellowstone' ts. metadata = list() #metadata$isRegularOrdered = TRUE # This arg not used any longer in this version metadata$whichDimIsTime = 1 # Which dim of the input refer to time for # 2D/3D inputs? Ignored for a single time # series input. metadata$startTime = c(1981,7,7) # Or startTime=1981.5137 # startTime=as.Date('1981-7-7') metadata$deltaTime = 1/24 # Half-monthly regular ts: 0.5/12=1/24 metadata$period = 1.0 # The period is 1 year: # freq x deltaTime = period # 24 x 1/24 = 1.0 metadata$omissionValue = NaN # By default, NaNs are ignored metadata$maxMissingRateAllowed = 0.7500 # If missingness is higher than .75, the ts # is skipped and not fitted metadata$deseasonalize = FALSE # Do not remove the global seasonal pattern # before fitting the beast model metadata$detrend = FALSE # Do not remove the global trend before # the fitting # prior is the ONLY true parameters of the beast model,used to specify the priors # in the Bayesian formulation prior = list() prior$seasonMinOrder = 1 #min harmonic order allowed to fit seasonal cmpnt prior$seasonMaxOrder = 5 #max harmonic order allowed to fit seasonal cmpnt prior$seasonMinKnotNum = 0 #min number of changepnts in seasonal cmpnt prior$seasonMaxKnotNum = 3 #max number of changepnts in seasonal cmpnt prior$seasonMinSepDist = 10 #min inter-chngpts separation for seasonal cmpnt prior$trendMinOrder = 0 #min polynomial order allowed to fit trend cmpnt prior$trendMaxOrder = 1 #max polynomial order allowed to fit trend cmpnt prior$trendMinKnotNum = 0 #min number of changepnts in trend cmpnt prior$trendMaxKnotNum = 15 #max number of changepnts in trend cmpnt prior$trendMinSepDist = 5 #min inter-chngpts separation for trend cmpnt prior$precValue = 10.0 #Initial value of the precision parameter (no # need to change it unless for precPrioType='const') prior$precPriorType = 'uniform' # Possible values: const, uniform, and componentwise # mcmc is NOT part of the beast model itself, but some parameters to configure the # MCMC inference. mcmc = list() mcmc$seed = 9543434# an arbitray seed for random number generator mcmc$samples = 3000 # samples collected per chain mcmc$thinningFactor = 3 # take every 3rd sample and discard others mcmc$burnin = 150 # discard the initial 150 samples per chain mcmc$chainNumber = 3 # number of chains mcmc$maxMoveStepSize = 4 # max random jump step when proposing new chngpts mcmc$trendResamplingOrderProb = 0.100 # prob of choosing to resample polynomial order mcmc$seasonResamplingOrderProb = 0.100 # prob of choosing to resample harmonic order mcmc$credIntervalAlphaLevel = 0.950 # the significance level for credible interval # extra is NOT part of the beast model itself, but some parameters to configure the # output and computation process extra = list() extra$dumpInputData = FALSE #If true, a copy of input time series is outputted extra$whichOutputDimIsTime = 1 #For 2D or 3D inputs, which dim of the output refers to # time? Ignored if the input is a single time series extra$computeCredible = FALSE #If true, compute CI: computing CI is time-intensive. extra$fastCIComputation = TRUE #If true, a faster way is used to get CI, but it is # still time-intensive. That is why the function beast() # is slow because it always compute CI. extra$computeSeasonOrder = FALSE #If true, dump the estimated harmonic order over time extra$computeTrendOrder = FALSE #If true, dump the estimated polynomial order over time extra$computeSeasonChngpt = TRUE #If true, get the most likely locations of s chgnpts extra$computeTrendChngpt = TRUE #If true, get the most likely locations of t chgnpts extra$computeSeasonAmp = FALSE #If true, get time-varying amplitude of seasonality extra$computeTrendSlope = FALSE #If true, get time-varying slope of trend extra$tallyPosNegSeasonJump= FALSE #If true, get those changpts with +/- jumps in season extra$tallyPosNegTrendJump = FALSE #If true, get those changpts with +/- jumps in trend extra$tallyIncDecTrendJump = FALSE #If true, get those changpts with increasing/ # decreasing trend slopes extra$printProgressBar = TRUE extra$printOptions = TRUE extra$quiet = FALSE # print warning messages, if any extra$consoleWidth = 0 # If 0, the console width is from the current console extra$numThreadsPerCPU = 2 # 'numThreadsPerCPU' and 'numParThreads' are used to extra$numParThreads = 0 # configure multithreading runs; they're used only if # Y has multiple time series (e.g.,stacked images) o = beast123(Yellowstone,metadata,prior,mcmc,extra, season='harmonic') plot(o) ## End(Not run) #------------------------------Example 4: Handle irregular time series-----------------# # Handle irregular time series: ohio is a data frame of a Landsat NDVI series observed # at unevely-spaced times ## Not run: data(ohio) str(ohio) metadata = list() metadata$time = ohio$time # Must supply individual times for irregular inputs metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$period = 1.0 o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### class(ohio$rdate) # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time = ohio$rdate # Must supply individual times for irregular inputs o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$Y # Another accepted time format for beast123 ohio$M ohio$M metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$year = ohio$Y metadata$time$month = ohio$M metadata$time$day = ohio$D o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$Y # Another accepted time format for beast123 ohio$doy metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$year = ohio$Y metadata$time$doy = ohio$doy o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$time # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time = ohio$time # Fractional year o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr1 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the time interval for aggregation metadata$time$datestr = ohio$datestr1 metadata$time$strfmt = '????yyyy?mm?dd' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr2 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply a desired time interval for aggregation metadata$time$datestr = ohio$datestr2 metadata$time$strfmt = '????yyyydoy????' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr3 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$datestr = ohio$datestr3 metadata$time$strfmt = 'Y,,M/D' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ## End(Not run) #------------------Example 4: Handle multiple time series (i.e., matrix input)-----------# # Handle multiple time series: 'simdata' is a 2D matrix of dim 300x3; it consits of 3 # time series of length 300 each. For this toy example, I decide to be lazy and use the same # time series for the three columns. ## Not run: data(simdata) # dim of simdata: 300 x 3 (time x num_of_time_series) dim(simdata) # the first dimenion refer to time (i.e, 300) metadata = list() metadata$whichDimIsTime = 1 # Which dim of the input refer to time for 2D inputs? # 300 is the ts length, so dim is set to '1' here. metadata$period = 24 # By default, we assume startTime=1 and deltaTime=1 extra=list() extra$whichOutputDimIsTime = 2 # Which dim of the output arrays refers to time? o=beast123(simdata, metadata,extra=extra) # Default values used for those missing parameters # The lists of arg parameters can also be directly provided inline within the command o=beast123( simdata, metadata=list(whichDimIsTime=1,period=24), extra=list(whichOutput=2) ) # The field names of the lists can be shortened as long as no ambiguitity is caused. o=beast123( simdata, metadata=list(whichDim=1,per=24), extra=list(whichOut=2) ) #------------------Example 4: Another run by transposing simdata--------------------------# simdata1=t(simdata) # dim of simdata1: 3 x 300 (num of ts x time ) metadata = list() metadata$whichDimIsTime = 2 # Which dim of the input refer to time for 2D inputs? # 300 is the ts length, so dim is set to '2' here. metadata$period = 24 # By default, we assume startTime=1 and deltaTime=1 o=beast123(simdata1, metadata) # Default values used for those missing parameters o=beast123( simdata1, metadata=list(whichDim=2, per=24) ) ## End(Not run) #------------------Example 5: Handle stacked time series images (e.g., 3d input)--------# # Handle 3D stacked images of irregular and unordered time-series: imagestack is a 3D # array of size 12x9x1066, each pixel being a time series of length 1066 ## Not run: data(imagestack) dim(imagestack$ndvi) # Dim: 12 x 9 X 1066 (row x col x time) imagestack$datestr # A character vector of 1066 date strings metadata = list() metadata$whichDimIsTime = 3 # Which dim of the input refer to time for 3D inputs? # 1066 is the ts length, so dim is set to '3' here. # In this example, this arg is not needed because # the time$datestr can also help to match and pick up # the right time dimesion of imagestack$ndvi. metadata$time$datestr = imagestack$datestr metadata$time$strfmt = 'LT05_018032_20080311.yyyy-mm-dd' metadata$deltaTime = 1/12 # Aggregate the irregular ts at a monthly interval:1/12 Yr metadata$period = 1.0 # The period is 1 year: deltaTime*freq=1/12*12=1.0 extra = list() extra$dumpInputData = TRUE # Get a copy of aggregated input ts extra$numThreadsPerCPU = 2 # Each cpu core will be assigned 2 threads extra$numParThreads = 0 # If 0, total_num_threads=numThreadsPerCPU*num_of_cpu_core # if >0, used to specify the total number of threads # Default values for missing parameters o=beast123(imagestack$ndvi, metadata=metadata,extra=extra) print(o,c(5,3)) # print the result for the pixel at Row 5 and Col 3 plot(o,c(5,3)) # plot the result for the pixel at Row 5 and Col 3 image(o$trend$ncp) # number of trend changepoints over space ## End(Not run) #---------Example 6: Handle stacked GeoTiff image files imported with the raster package------# # Handle 3D stacked images of irregular time-series : 'ndvi.zip' is a zip file of # 437 NDIV tiff image files, each having a dim of 12 x 9. # Code availlable at https://github.com/zhaokg/Rbeast/blob/master/R/beast123_raster_example.txt#--------------------------------NOTE-------------------------------------------------# # beast123() is an all-inclusive function that duplicates the functionalities of beast # and beast.irreg. It can handle a single, multiple, or 3D of stacked time series, being # either regular or irregular. It allows for customization through four LIST arguments: # metadata -- additional info about the input Y # prior -- prior parameters for the beast model # mcmc -- MCMC simulation setting # extra -- misc parameters turning on/off outputs and setting up parallel computations # # Despite being essentially the same as beast and beast.irreg, beast123 is provided mainly # to support concurrent handling of multiple time series (e.g., stacked satellite images) # via parallel computing: When processing stacked raster layers, DO NOT iterate pixel by pixel # using beast() or beast.irreg() via an external parallel caller (e.g., doParallel or foreach). # Instread, please use beast123(), which supports mulithreading internally. #------------------------------Example 1: one time series with seasonalty-------------# # Yellowstone is a half-monthly time series of 774 NDVI measurmments at a Yellowstone # site starting from July 1-15,1981(i.e., start=c(1981,7,7). It has 24 data points per # year (freq=24). library(Rbeast) data(Yellowstone) plot(Yellowstone) # Below, the four option args are missing, so defalut values will be used, with some # warning messages given to altert this. By default, the input Y is assumed to be regular # with a seasonal component. The default arg values used will be printed out and they can # serve as a template to customize the parameters. o = beast123(Yellowstone) plot(o) #------------------------------Example 2: a trend-only time series-------------------# # Nile is an annual river flow time series (i.e., no periodic variation). So, season # is set to 'none' to indicate trend-only analysis. Default values are used for other # missing options. Unlike the beast() function, beast123 does NOT use the time attributes # of a 'ts' object. For example, Nile is treated as a pure data number; its (start=1871, # end=1970, freq=1) attributes are ignored. The default times 1:length(Nile) are used # instead. The true time info need to be specified by the 'metadata' parameter, as shown # in the next example. o = beast123(Nile,season='none') plot(o) #------------------------------Example 3: call via the full API interface-----------# # Specify metadata, prior, mcmc, and extra explicitly. Only 'prior' is the true statistical # model parameters of BEAST; the other three are just options to configure the input/ouput # or the computation process. ## Not run: # metadata is NOT part of BEAST itself, but some extra info to describe the input # time series Y. Below, the input Y is the 'Yellowstone' ts. metadata = list() #metadata$isRegularOrdered = TRUE # This arg not used any longer in this version metadata$whichDimIsTime = 1 # Which dim of the input refer to time for # 2D/3D inputs? Ignored for a single time # series input. metadata$startTime = c(1981,7,7) # Or startTime=1981.5137 # startTime=as.Date('1981-7-7') metadata$deltaTime = 1/24 # Half-monthly regular ts: 0.5/12=1/24 metadata$period = 1.0 # The period is 1 year: # freq x deltaTime = period # 24 x 1/24 = 1.0 metadata$omissionValue = NaN # By default, NaNs are ignored metadata$maxMissingRateAllowed = 0.7500 # If missingness is higher than .75, the ts # is skipped and not fitted metadata$deseasonalize = FALSE # Do not remove the global seasonal pattern # before fitting the beast model metadata$detrend = FALSE # Do not remove the global trend before # the fitting # prior is the ONLY true parameters of the beast model,used to specify the priors # in the Bayesian formulation prior = list() prior$seasonMinOrder = 1 #min harmonic order allowed to fit seasonal cmpnt prior$seasonMaxOrder = 5 #max harmonic order allowed to fit seasonal cmpnt prior$seasonMinKnotNum = 0 #min number of changepnts in seasonal cmpnt prior$seasonMaxKnotNum = 3 #max number of changepnts in seasonal cmpnt prior$seasonMinSepDist = 10 #min inter-chngpts separation for seasonal cmpnt prior$trendMinOrder = 0 #min polynomial order allowed to fit trend cmpnt prior$trendMaxOrder = 1 #max polynomial order allowed to fit trend cmpnt prior$trendMinKnotNum = 0 #min number of changepnts in trend cmpnt prior$trendMaxKnotNum = 15 #max number of changepnts in trend cmpnt prior$trendMinSepDist = 5 #min inter-chngpts separation for trend cmpnt prior$precValue = 10.0 #Initial value of the precision parameter (no # need to change it unless for precPrioType='const') prior$precPriorType = 'uniform' # Possible values: const, uniform, and componentwise # mcmc is NOT part of the beast model itself, but some parameters to configure the # MCMC inference. mcmc = list() mcmc$seed = 9543434# an arbitray seed for random number generator mcmc$samples = 3000 # samples collected per chain mcmc$thinningFactor = 3 # take every 3rd sample and discard others mcmc$burnin = 150 # discard the initial 150 samples per chain mcmc$chainNumber = 3 # number of chains mcmc$maxMoveStepSize = 4 # max random jump step when proposing new chngpts mcmc$trendResamplingOrderProb = 0.100 # prob of choosing to resample polynomial order mcmc$seasonResamplingOrderProb = 0.100 # prob of choosing to resample harmonic order mcmc$credIntervalAlphaLevel = 0.950 # the significance level for credible interval # extra is NOT part of the beast model itself, but some parameters to configure the # output and computation process extra = list() extra$dumpInputData = FALSE #If true, a copy of input time series is outputted extra$whichOutputDimIsTime = 1 #For 2D or 3D inputs, which dim of the output refers to # time? Ignored if the input is a single time series extra$computeCredible = FALSE #If true, compute CI: computing CI is time-intensive. extra$fastCIComputation = TRUE #If true, a faster way is used to get CI, but it is # still time-intensive. That is why the function beast() # is slow because it always compute CI. extra$computeSeasonOrder = FALSE #If true, dump the estimated harmonic order over time extra$computeTrendOrder = FALSE #If true, dump the estimated polynomial order over time extra$computeSeasonChngpt = TRUE #If true, get the most likely locations of s chgnpts extra$computeTrendChngpt = TRUE #If true, get the most likely locations of t chgnpts extra$computeSeasonAmp = FALSE #If true, get time-varying amplitude of seasonality extra$computeTrendSlope = FALSE #If true, get time-varying slope of trend extra$tallyPosNegSeasonJump= FALSE #If true, get those changpts with +/- jumps in season extra$tallyPosNegTrendJump = FALSE #If true, get those changpts with +/- jumps in trend extra$tallyIncDecTrendJump = FALSE #If true, get those changpts with increasing/ # decreasing trend slopes extra$printProgressBar = TRUE extra$printOptions = TRUE extra$quiet = FALSE # print warning messages, if any extra$consoleWidth = 0 # If 0, the console width is from the current console extra$numThreadsPerCPU = 2 # 'numThreadsPerCPU' and 'numParThreads' are used to extra$numParThreads = 0 # configure multithreading runs; they're used only if # Y has multiple time series (e.g.,stacked images) o = beast123(Yellowstone,metadata,prior,mcmc,extra, season='harmonic') plot(o) ## End(Not run) #------------------------------Example 4: Handle irregular time series-----------------# # Handle irregular time series: ohio is a data frame of a Landsat NDVI series observed # at unevely-spaced times ## Not run: data(ohio) str(ohio) metadata = list() metadata$time = ohio$time # Must supply individual times for irregular inputs metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$period = 1.0 o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### class(ohio$rdate) # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time = ohio$rdate # Must supply individual times for irregular inputs o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$Y # Another accepted time format for beast123 ohio$M ohio$M metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$year = ohio$Y metadata$time$month = ohio$M metadata$time$day = ohio$D o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$Y # Another accepted time format for beast123 ohio$doy metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$year = ohio$Y metadata$time$doy = ohio$doy o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$time # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time = ohio$time # Fractional year o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr1 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the time interval for aggregation metadata$time$datestr = ohio$datestr1 metadata$time$strfmt = '????yyyy?mm?dd' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr2 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply a desired time interval for aggregation metadata$time$datestr = ohio$datestr2 metadata$time$strfmt = '????yyyydoy????' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ##################################################################################### ohio$datestr3 # Another accepted time format for beast123 metadata = list() metadata$deltaTime = 1/12 # Must supply the desired time interval for aggregation metadata$time$datestr = ohio$datestr3 metadata$time$strfmt = 'Y,,M/D' o=beast123(ohio$ndvi, metadata) # Default values used for those missing parameters ## End(Not run) #------------------Example 4: Handle multiple time series (i.e., matrix input)-----------# # Handle multiple time series: 'simdata' is a 2D matrix of dim 300x3; it consits of 3 # time series of length 300 each. For this toy example, I decide to be lazy and use the same # time series for the three columns. ## Not run: data(simdata) # dim of simdata: 300 x 3 (time x num_of_time_series) dim(simdata) # the first dimenion refer to time (i.e, 300) metadata = list() metadata$whichDimIsTime = 1 # Which dim of the input refer to time for 2D inputs? # 300 is the ts length, so dim is set to '1' here. metadata$period = 24 # By default, we assume startTime=1 and deltaTime=1 extra=list() extra$whichOutputDimIsTime = 2 # Which dim of the output arrays refers to time? o=beast123(simdata, metadata,extra=extra) # Default values used for those missing parameters # The lists of arg parameters can also be directly provided inline within the command o=beast123( simdata, metadata=list(whichDimIsTime=1,period=24), extra=list(whichOutput=2) ) # The field names of the lists can be shortened as long as no ambiguitity is caused. o=beast123( simdata, metadata=list(whichDim=1,per=24), extra=list(whichOut=2) ) #------------------Example 4: Another run by transposing simdata--------------------------# simdata1=t(simdata) # dim of simdata1: 3 x 300 (num of ts x time ) metadata = list() metadata$whichDimIsTime = 2 # Which dim of the input refer to time for 2D inputs? # 300 is the ts length, so dim is set to '2' here. metadata$period = 24 # By default, we assume startTime=1 and deltaTime=1 o=beast123(simdata1, metadata) # Default values used for those missing parameters o=beast123( simdata1, metadata=list(whichDim=2, per=24) ) ## End(Not run) #------------------Example 5: Handle stacked time series images (e.g., 3d input)--------# # Handle 3D stacked images of irregular and unordered time-series: imagestack is a 3D # array of size 12x9x1066, each pixel being a time series of length 1066 ## Not run: data(imagestack) dim(imagestack$ndvi) # Dim: 12 x 9 X 1066 (row x col x time) imagestack$datestr # A character vector of 1066 date strings metadata = list() metadata$whichDimIsTime = 3 # Which dim of the input refer to time for 3D inputs? # 1066 is the ts length, so dim is set to '3' here. # In this example, this arg is not needed because # the time$datestr can also help to match and pick up # the right time dimesion of imagestack$ndvi. metadata$time$datestr = imagestack$datestr metadata$time$strfmt = 'LT05_018032_20080311.yyyy-mm-dd' metadata$deltaTime = 1/12 # Aggregate the irregular ts at a monthly interval:1/12 Yr metadata$period = 1.0 # The period is 1 year: deltaTime*freq=1/12*12=1.0 extra = list() extra$dumpInputData = TRUE # Get a copy of aggregated input ts extra$numThreadsPerCPU = 2 # Each cpu core will be assigned 2 threads extra$numParThreads = 0 # If 0, total_num_threads=numThreadsPerCPU*num_of_cpu_core # if >0, used to specify the total number of threads # Default values for missing parameters o=beast123(imagestack$ndvi, metadata=metadata,extra=extra) print(o,c(5,3)) # print the result for the pixel at Row 5 and Col 3 plot(o,c(5,3)) # plot the result for the pixel at Row 5 and Col 3 image(o$trend$ncp) # number of trend changepoints over space ## End(Not run) #---------Example 6: Handle stacked GeoTiff image files imported with the raster package------# # Handle 3D stacked images of irregular time-series : 'ndvi.zip' is a zip file of # 437 NDIV tiff image files, each having a dim of 12 x 9. # Code availlable at https://github.com/zhaokg/Rbeast/blob/master/R/beast123_raster_example.txt
DNA copy number alteration data in array-based CGH data for Chromesome 11
Description
CNAchrom11 is a vector of the log2 intensity ratios for cell line GM03576 for Chromosome 11, obtained from Snijders et al. (2001).
Usage
data(CNAchrom11)data(CNAchrom11)
Source
Snijders et al. (2001), Assembly of microarrays for genome-wide measurement of DNA copy number, Nature Genetics, 29, 263-264 (http://www.nature.com/ng/journal/v29/n3/full/ng754.html).
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(CNAchrom11) o = beast(CNAchrom11, season='none') # no periodic component plot(o)library(Rbeast) data(CNAchrom11) o = beast(CNAchrom11, season='none') # no periodic component plot(o)
Daily confirmed COVID19 cases and deaths in the world
Description
covid19 is a data frame consisting of daily confirmed COVID19 cases and deaths in the world from Jan 22, 2020 to Dec 16, 2021.
Usage
data(covid19)data(covid19)
Source
https://ourworldindata.org/grapher/daily-covid-cases-deaths?country=~OWID_WRL (last accessed on Dec 16, 2021)
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(covid19) plot(covid19$date, covid19$newcases, type='l') ## Not run: # Apply a square root-transformation newcases = sqrt( covid19$newcases ) # This time series varies periodically every 7 days. 7 days can't be precisely # represented in the unit of year bcz some years has 365 days and others has 366. # BEAST can hanlde this in two ways. #(1) Use the date number as the time unit--the num of days lapsed since 1970-01-01. datenum = as.numeric(covid19$date) o = beast(newcases, start=min(datenum), deltat=1, period=7) o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date. plot(o) #(2) Use strings to explicitly specify deltat and period with a unit. startdate = covid19$date[1] o = beast(newcases, start=startdate, deltat='1day', period='7days') plot(o) ## End(Not run)library(Rbeast) data(covid19) plot(covid19$date, covid19$newcases, type='l') ## Not run: # Apply a square root-transformation newcases = sqrt( covid19$newcases ) # This time series varies periodically every 7 days. 7 days can't be precisely # represented in the unit of year bcz some years has 365 days and others has 366. # BEAST can hanlde this in two ways. #(1) Use the date number as the time unit--the num of days lapsed since 1970-01-01. datenum = as.numeric(covid19$date) o = beast(newcases, start=min(datenum), deltat=1, period=7) o$time = as.Date(o$time, origin='1970-01-01') # Convert from integers to Date. plot(o) #(2) Use strings to explicitly specify deltat and period with a unit. startdate = covid19$date[1] o = beast(newcases, start=startdate, deltat='1day', period='7days') plot(o) ## End(Not run)
Landsat reflectance and NDVI time series from Google Earth Engine
Description
Get Landsat reflectance and NDVI time series from Google Earth Engine given longitude and latitude
Usage
geeLandsat(lon=NA, lat=NA, radius=100, stat='mean',timeout=700)geeLandsat(lon=NA, lat=NA, radius=100, stat='mean',timeout=700)
Arguments
lon |
numeric within [-180,180] |
lat |
numeric within [-90, 90] |
radius |
a positive number ( <=500 meters ); the radius of a buffer around the given latitude and longitude for aggregation. If |
stat |
character; if |
timeout |
integer; the seconds elapsed to wait for connection timeout. See the note for an explanation. |
Value
a data.frame object consisting of dates, sensor type, reflectances, and NDVI for the requested location. It contains only valid and clear-sky values as obtained by referring to the standard clouds flags.
Note
As a poor man's scheme to interact with Google Earth Engine, geeLandsat should be used only for occasional retrieval of Landsat time series at a few sites, NOT for batch downloading for thousands of sites in a R loop. This procedure is provided to get example time series for testing BEAST. Behind the scene, this function calls to a free Python-based server using my own GEE credential. Normally it takes several seconds to retrieve one time series, but as a free cloud service, the Python server only offers 100 seconds of free CPU time per day, with throttling applied. So it may take up to a few mins to get a time series on your end. It may fail due to connection timeout; if so, give it a few tries. If you need to retrieve data for thousands or millions of sites, please contact the author.
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, minesweeper, tetris
Examples
library(Rbeast) ## Not run: df = geeLandsat(lon=-80.983877,lat= 40.476882) #if it fails, try a few more times before giving up print(df) ## End(Not run)library(Rbeast) ## Not run: df = geeLandsat(lon=-80.983877,lat= 40.476882) #if it fails, try a few more times before giving up print(df) ## End(Not run)
A monthly Google Trend time series of the US search interest in the word "beach"
Description
googletrend_beach is a ts object comprising monthly search interest in "beach" from the United States, as reported from Google Trends. Sudden changes in the search trend are attributed to extreme weather events or the covid19 outbreak
Usage
data(googletrend_beach)data(googletrend_beach)
Source
https://trends.google.com/trends/explore?date=all&geo=US&q=beach
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(googletrend_beach) # A monthly ts starting from Jan 2004 o = beast(googletrend_beach ) plot(o)library(Rbeast) data(googletrend_beach) # A monthly ts starting from Jan 2004 o = beast(googletrend_beach ) plot(o)
Decades of Landsat NDVI time series over a small area in Ohio
Description
imagestack is a LIST containing Landsat-derived NDVI image chips at an Ohio site
Usage
data(imagestack)data(imagestack)
Source
Landsat images courtesy of the U.S. Geological Survey
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
data(imagestack) imagestack$datestr # A string vector containg the observation dates of individual ndvi images ## Not run: imagestack$ndvi # NDVI images collected over the past several deccades ## End(Not run) plot(imagestack$ndvi[3,4,],type='l') # Plot the raw data at a pixeldata(imagestack) imagestack$datestr # A string vector containg the observation dates of individual ndvi images ## Not run: imagestack$ndvi # NDVI images collected over the past several deccades ## End(Not run) plot(imagestack$ndvi[3,4,],type='l') # Plot the raw data at a pixel
The Minesweeper game in R
Description
A poor man's implementation of the minesweeper game in R. Yes, you are right: it has nothing to do with time series decomposition, changepoint detection, and time series segmentation. Its only remote connection to Rbeast is that this is a practice script I wrote to learn R graphics for implementing Rbeast.
Usage
minesweeper(height=15, width=12, prob=0.1)minesweeper(height=15, width=12, prob=0.1)
Arguments
height |
integer; number of rows of the mine grid along the vertical direction. |
width |
integer; number of columns of the mine grid along the horizontal direction. |
prob |
numeric; a fraction between 0 and 1 to specify the probability of mine occurrence in the mine grid. |
Value
Instructions:
-
LEFT-click to clear a spot.
-
RIGHT-click to flag a spot.
-
MIDDLE-click(wheel) a cleared and numbered spot to open neighbor spots, if flagged correctly.
-
Click Restart for a new game
Note
An interactive graphics window is needed to run this function correctly. So it won't run in RStudio's plot pane. The function will use the x11() or x11(type='Xlib') graphic device to open a pop-up window.
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, tetris, geeLandsat
Examples
library(Rbeast) ## Not run: minesweeper() # A mine field of size 20x25 with rougly a 15 minesweeper(20,25,0.15) ## End(Not run)library(Rbeast) ## Not run: minesweeper() # A mine field of size 20x25 with rougly a 15 minesweeper(20,25,0.15) ## End(Not run)
An irregular Landsat NDVI time series at an Ohio site
Description
ohio is a data.frame object comprising decades of Landsat-observed surface reflectances and NDVI at an Ohio site
Usage
data(ohio)data(ohio)
Source
Landsat images courtesy of the U.S. Geological Survey
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(ohio) # Landsat surface references and NDVI at a single pixel observed over time str(ohio) ## Not run: # ohio$ndvi is a single irregular time series y = ohio$ndvi o = beast.irreg(y, time=ohio$time,deltat=1/12) plot(o) print(o) # ohio also contains irregular time series of individual spectral bands # Below, run the multivariate version of the BEAST algorithm to decompose # the 5 time series and detect common changepoints altogether y = list(ohio$blue, ohio$green, ohio$red, ohio$nir, ohio$swir1); o = beast.irreg(y, time=ohio$time,deltat=1/12, freq=12) plot(o) print(o) ## End(Not run)library(Rbeast) data(ohio) # Landsat surface references and NDVI at a single pixel observed over time str(ohio) ## Not run: # ohio$ndvi is a single irregular time series y = ohio$ndvi o = beast.irreg(y, time=ohio$time,deltat=1/12) plot(o) print(o) # ohio also contains irregular time series of individual spectral bands # Below, run the multivariate version of the BEAST algorithm to decompose # the 5 time series and detect common changepoints altogether y = list(ohio$blue, ohio$green, ohio$red, ohio$nir, ohio$swir1); o = beast.irreg(y, time=ohio$time,deltat=1/12, freq=12) plot(o) print(o) ## End(Not run)
Bayesian changepoint detection and time series decomposition
Description
Plot the result obtained from the beast function.
Usage
## S3 method for class 'beast' plot( x, index = 1, vars = c('y','s','scp','sorder','t','tcp','torder','slpsgn','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'), ... )## S3 method for class 'beast' plot( x, index = 1, vars = c('y','s','scp','sorder','t','tcp','torder','slpsgn','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'), ... )
Arguments
x |
a "beast" object returned by |
index |
an integer (default to 1 ) or a vector of two integers to specify the index of the time series to plot if |
vars |
a vector of strings indicating the elements or variables of |
relative.heights |
a numeric vector of the same length as that of |
col |
a string vector of the same length as that of |
main |
a string; the main title. |
xlab |
a string: the x axis title. |
ylab |
a string vector of the same length as that of |
cex.main |
cex for the main title |
cex.lab |
cex for the axis title |
interactive |
a bool scalar. If TRUE, an interactive GUI is used for examining individual elements of |
ncpStat |
character. A string to specify which statistic is used for the Number of ChangePoint (ncp). Five values are possible: 'mean', 'mode', 'median','pct90', and 'max'; the default is 'median'. Individual models sampled by BEAST has a varying dimension (e.g., number of changepoints or knots). For example, if mcmc$samples=10, the numbers of changepoints for the 10 sampled models are assumed to be c(0, 2, 4, 1, 1, 2, 7, 6, 6, 1). The mean ncp will be 3.1 (rounded to 3), the median is 2.5 (2), the mode is 1, and the maximum is 7. The 'max' option plots all the changepoints recorded in |
... |
additional parameters to be implemented. |
Value
This function creates various plots to demonstrate the results of a beast decomposition. .
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, plot.beast,minesweeper, tetris, geeLandsat
Examples
library(Rbeast) data(simdata) ## Not run: result=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(result,1) plot(result,2) ## End(Not run)library(Rbeast) data(simdata) ## Not run: result=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(result,1) plot(result,2) ## End(Not run)
Bayesian changepoint detection and time series decomposition
Description
Summarize and print the results obtained from the BEAST time series decomposition and segmentation.
Usage
## S3 method for class 'beast' print( x, index = 1, ... )## S3 method for class 'beast' print( x, index = 1, ... )
Arguments
x |
a "beast" object returned by |
index |
an integer (default to 1 ) or a vector of two integers to specify the index of the time series to print if |
... |
additional parameters to be implemented. |
Value
Print a summary of changepoints detected for the seasonal or trend component.
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, minesweeper, tetris, geeLandsat
Examples
library(Rbeast) data(simdata) ## Not run: #out=beast123(simdata) #Error: whichDimIsTime has to be specified to # tell which dim of simdata refers to time. # See below. out=beast123(simdata, metadata=list(whichDimIsTime=1)) print(out, 1) print(out, 2) ## End(Not run)library(Rbeast) data(simdata) ## Not run: #out=beast123(simdata) #Error: whichDimIsTime has to be specified to # tell which dim of simdata refers to time. # See below. out=beast123(simdata, metadata=list(whichDimIsTime=1)) print(out, 1) print(out, 2) ## End(Not run)
Simulated time series to test BEAST
Description
simdata is a 300 x 3 matrix, consisting three time series of length 300. Currently, the three time series are the same. It is used to illustrate BEAST can handle multiple time series at a single function call.
of BEAST.
Usage
data(simdata)data(simdata)
Source
Rbeast v0.9.2
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(simdata) plot(simdata[,1],type='l') ## Not run: #out=beast123(simdata) # Error: whichDimIsTime has to be specified. See below out=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(out,1) plot(out,2) plot(out,3) ## End(Not run)library(Rbeast) data(simdata) plot(simdata[,1],type='l') ## Not run: #out=beast123(simdata) # Error: whichDimIsTime has to be specified. See below out=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(out,1) plot(out,2) plot(out,3) ## End(Not run)
The Tetris game in R
Description
A poor man's implementation of the Tetris game in R. Yes, you are right again: it has nothing to do with time series decomposition, changepoint detection, and time series segmentation. Its only remote connection to Rbeast is that this is a practice script I wrote to learn R graphics for implementing Rbeast.
Usage
tetris(height=25, width=14, speed=0.6)tetris(height=25, width=14, speed=0.6)
Arguments
height |
integer; number of rows of the mine grid along the vertical direction. |
width |
integer; number of columns of the mine grid along the horizontal direction. |
speed |
numeric; a time interval between 0.05 and 2 seconds, specifying how fast the tetriminos moves down. The smaller, the faster. |
Value
Instructions:
-
Left arrow to move left.
-
Right arrow to move right.
-
Up arrow to rotate.
-
Down arrow to speed up.
-
Space key to sink to the bottom.
Note
This function works only under the Windows OS not Linux or Mac. An interactive graphics window is needed to run this function correctly. So it won't run in RStudio's plot pane. The function will use the x11() or x11(type='Xlib') graphic device to open a pop-up window.
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, minesweeper, geeLandsat
Examples
library(Rbeast) ## Not run: tetris() # A field of size 20x25 with blocks moving down every 0.1 sec. tetris(20,25,0.1) ## End(Not run)library(Rbeast) ## Not run: tetris() # A field of size 20x25 with blocks moving down every 0.1 sec. tetris(20,25,0.1) ## End(Not run)
Bayesian changepoint detection and time series decomposition
Description
Extract the result of a single time series from an object of class beast
Usage
tsextract( x, index = 1 )tsextract( x, index = 1 )
Arguments
x |
a "beast" object returned by |
index |
an integer (default to 1 ) or a vector of two integers to specify the index of the time series to extract if |
Value
A LIST object of the result for the chosen time series, which contains the same field as x.
Note
Use this function only to manually and interactively examine individual times series. If the purpose is to loop through x, the use of direct indexing is much faster. For example, if x is a beast object for a 300x200x1000 3D array (row x col x time), use x$trend$Y[20,40,] to get the fitted trend at the pixel of row 20 and col 40.
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
See Also
beast, beast.irreg, beast123, minesweeper, tetris, geeLandsat
Examples
library(Rbeast) data(simdata) # handle only the 1st ts out=beast(simdata[,1]) ## Not run: # handle all the ts out=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(out,1) plot(out,2) ## End(Not run)library(Rbeast) data(simdata) # handle only the 1st ts out=beast(simdata[,1]) ## Not run: # handle all the ts out=beast123(simdata, metadata=list(whichDimIsTime=1)) plot(out,1) plot(out,2) ## End(Not run)
30 years' AVHRR NDVI data at a Yellostone site
Description
Yellowstone is a vector comprising 30 years' AVHRR NDVI data at a Yellostone site
Usage
data(Yellowstone)data(Yellowstone)
Source
Rbeast v0.9.2
References
-
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, 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.102-119 (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 fine-scale human disturbances in a working landscape with Landsat time series on Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 176, pp.250-261(a beast application paper).
Examples
library(Rbeast) data(Yellowstone) plot(Yellowstone,type='l') result=beast(Yellowstone) plot(result)library(Rbeast) data(Yellowstone) plot(Yellowstone,type='l') result=beast(Yellowstone) plot(result)