#--------------------------------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.
# \donttest{
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.
# \donttest{
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.
if (FALSE) {
# 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 # Regular input
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$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)
}
#------------------------------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
# \donttest{
data(ohio)
str(ohio)
metadata = list()
metadata$isRegularOrdered = FALSE # The input data is irregular
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$isRegularOrdered = FALSE # The input data is irregular
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$isRegularOrdered = FALSE # The input data is irregular
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$isRegularOrdered = FALSE # Irregular input
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$isRegularOrdered = FALSE # Irregular input
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$isRegularOrdered = FALSE # Irregular input
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$isRegularOrdered = FALSE # Irregular input
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$isRegularOrdered = FALSE # Irregular input
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
# }
#------------------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.
if (FALSE) {
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) )
}
#------------------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
if (FALSE) {
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$isRegularOrdered = FALSE # 'imagestack$ndvi' is an IRREGULAR input
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.
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
}
#---------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
Run the code above in your browser using DataLab