# NOT RUN {
library(Rbeast)
# A MODIS time series of NDVI for a forest plot in Ohio. It has 23 samples
# per year (i.e., period=23). Note that the input time series to "beast" must
# be spaced/observed at regular time intervals, with missing values indicated
# by NAs, NaNs, or a custimized value (see Example 3). Iregular-sampled time
# series data need to be first aggegrated at a regular time interval of your
# choice before running beast; if not, beast may give meaningless results or,
# even worse, terminate abnormally.
data(modis_ohio)
plot(modis_ohio)
# }
# NOT RUN {
#--------------------------------Example 1--------------------------------#
# No "option" argument supplied below so default parameters are used. The period
# (i.e., 23) will be estimated via auto-correlation. Letting the program compute the
# period of a cyclic time series is not always reliable, so it is always suggested
# to directly supply the period as in Example 2 and Example 3.
out=beast(modis_ohio)
plot(out)
#*****************************End of Example 1****************************#
#--------------------------------Example 2--------------------------------#
# "option" is set to 23, specicfying the period of modis_ohio as 23
out=beast(modis_ohio,23)
plot(out)
plot(out$s) #The same as plot(out$s[,1]): plot the seasonal curve
plot(out$sProb) #Plot the probability of observing seasonal changepoints
plot(out$t) #The same as plot(out$t[,1]): plot the trend
plot(out$sProb) #Plot the probability of observing trend changepoints
# }
# NOT RUN {
#*****************************End of Example 2****************************#
#--------------------------------Example 3--------------------------------#
# Specify the option parameters explicilty
opt=list() #Create an empty list to append individual model parameters
opt$period=23 #Period of the cyclic/seasonal component of the modis time series
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 seperation time btw neighboring season change-pts(mustbe >=0)
opt$minSepDist_Trend=20 #Min seperation time btw neighboring trend change-pts(must be >=0)
opt$maxKnotNum_Season=4 #Max number of season changepoints allowed
opt$maxKnotNum_Trend=10 #Max number of trend changepoints allowed
opt$omittedValue=-999 #A customized value to indicate bad/missing values in the time
#series, in additon to those NA or NaN values.
opt$printToScreen=1 #If set to 1, display some progress status while running
opt$printCharLen=150 #The length of chars in each status line when printToScreen=1
# 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 (e.g., samples taken every 3 iterations)
opt$burnin=500 #Number of burn-in samples discarded at the start of each chain
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.
opt$computeCredible=0 #If set to 1, compute 95% credible intervals: The results will be
#saved as sCI, tCI, and bCI in the output variable.
opt$fastCIComputation=0 #If set to 1, employ a fast algorithm to compute credible intervals
opt$computeSlopeSign=1 #If set to 1, compute the probability of having a postive slope in
#the estimated trend. The result will be saved as bsign in the output
#variable.
opt$computeHarmonicOrder=1 #If set to 1, compute the mean harmonic order of the fitted
#seasonal component. The result will be saved as "horder" in
#the output variable.
opt$computeTrendOrder=1 #If set to 1, compute the mean polynomial order of the fitted
#trend component. The result will be saved as "torder" in
#the output variable.
#opt$outputToDisk=0 #(if set to 1, results will be written to files in a folder)
#opt$outputFolder ='c:/out'#Specify the output folder when outputToDisk=1
#opt$lengthPerTimeSeries_infile=300#the time series length if input data come from a binary file
# }
# NOT RUN {
# Use "opt" defined above in the beast function. Note that to run beast(), not all the individual
# parameters in option need to be explicitly specified. If an parameter is not given in option, its
# default value will be used.
out=beast(modis_ohio, opt)
plot(out)
# }
# NOT RUN {
#*****************************End of Example 3****************************#
# }
# NOT RUN {
#--------------------------------Example 4--------------------------------#
# 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.
beast(modis_ohio, 23, demoGUI=TRUE)
#*****************************End of Example 4****************************#
#--------------------------------Example 5--------------------------------#
# 'simdata' is a 300x3 matrix, consisting of three time series
data(simdata)
# Plot individual time series. As a toy example, all the three time series
# are the same.
plot(simdata[,1])
plot(simdata[,2])
# Below, the option is defined in the command line as a temporary list.
out=beast( simdata, list(period=24, chainNumber=3, sample=1000, burnin=200) )
# "out" contains results for the three time series. Plot the result for the second one
plot(out,2)
#*****************************End of Example 5****************************#
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab