Skip to content

Instantly share code, notes, and snippets.

@danhammer
Created January 28, 2012 00:05
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save danhammer/1691694 to your computer and use it in GitHub Desktop.
Save danhammer/1691694 to your computer and use it in GitHub Desktop.
trend analysis in python and r
def convert_ts(time_series, start_year=2000, start_pd=4, freq=23):
"""
Convert a numpy time-series into an rpy2 object, which, in turn,
is a 'ts' object in R. The 'ts' object is more wholly specified
if a start date is provided. For our applications at 16-day
intervals for NDVI, this start date is April 2000, with a
frequency of 23 observations each year.
input: numpy time-series
output: rpy2 ts object
"""
r_ts = robjects.r.ts(robjects.FloatVector(time_series),
start=robjects.r.c(start_year,start_pd), frequency=freq)
return r_ts
def seasonal_decomposition_R(time_series):
"""
Seasonal decomposition of time-series by Loess smoothing ("stl"
function in R package "stats"). We remove the seasonal component
using the window option "periodic" whereby smoothing is
effectively replaced by taking the mean (see documentation for
"stl"). This is the method used for the BFAST algorithm.
input: rpy2 object ('ts' object in R)
output: numpy array
"""
stats = importr('stats')
robjects.globalenv["rts"] = convert_ts(time_series)
ts = robjects.r('c(rts - stl(rts, "periodic")$time.series[, "seasonal"])')
return np.array(ts)
def break_statistic(time_series):
"""
Test whether there exists a break in the time_series.
input: numpy array of time-series values
output: float value to test against, say, p=0.05
"""
bfast = importr('bfast')
robjects.globalenv["rts"] = convert_ts(time_series)
robjects.r('sct <- sctest(efp(rts ~ time(rts), h=0.05, type="OLS-MOSUM"))')
return robjects.r('sct$statistic')[0]
@BastosA
Copy link

BastosA commented Apr 30, 2017

Hello

thanks!
I am running python 2.7 on anaconda. I added the bfast package to the R libraries in anaconda, but when I load it the script crashes. Any idea what might be happening?
I tried loading bfast from another folder, but it doesn't seem to be recognized...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment