Start using git to manage the source code

.Rbuildignore

+^.*\.Rproj$
+^\.Rproj\.user$
+R/RcppExports.R

.gitignore

+.Rproj.user
+.Rhistory
+.RData
+Langevin.Rproj
+src/*.o
+src/*.so
+src/*.dll

ChangeLog

+2016-01-27  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * OpenMP: Use recommended method to set the number of used threads
+
+2015-11-02  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add missing importFroms
+    * Tag as version 1.1.1
+
+2015-11-02  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add citation file
+    * Update vignette
+    * Tag as version 1.1
+
+2015-10-21  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add plot function for class 'Langevin' (only for 1D)
+
+2015-10-19  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add print function for class 'Langevin'
+
+2015-10-16  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add summary function for class 'Langevin'
+
+2015-09-07  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin{1,2}D: read sampling frequency from time-series object
+        * Langevin2D: time-series need to be arranged as columns now!
+
+2015-09-03  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries{1,2}D: output is a time-series object now
+        * timeseries2D: time-series are the columns of the returned object now!
+
+2015-08-26  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Use title case for package title
+    * Tag as version 1.0.3
+
+2015-08-26  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Don't install vignette source files in inst/doc
+    * Tag as version 1.0.2
+
+2015-08-25  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add paper submitted to JSS as vignette
+    * Tag as version 1.0.1
+
+2015-08-25  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Update package description
+    * Tag as version 1.0
+
+2015-07-21  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin1D: calculate error of the Diffusion
+
+2015-07-21  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add description of package
+    * Langevin1D: move part of the documentation to package description
+
+2015-07-16  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries2D: major cleanup
+        * rewrite to use matrices as input for the coefficients
+        * use cubic drift and quadratic diffusion polynomial
+        * output is a matrix with two rows now
+    * timeseries1D: use cubic drift and quadratic diffusion polynomial
+
+2015-07-15  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Rename timeseries -> timeseries1D to be consistent with Langevin{1,2}D
+    * timeseries{1,2}D: remove argument seed:
+        * One should use set.seed(seed) before calling timeseries{1,2}D instead
+
+2015-02-16  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Use 'Rcpp Attributes' to generate the glue between C++ and R.
+
+2015-01-13  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add Pedro G. Lind to package authors
+
+2015-01-13  Pedro G. Lind <pedro.g.lind@forwind.de>
+    * Add function timeseries2D: Generate a two-dimensional Langevin process
+
+2014-09-03  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin2D: Do the linear regression right
+    * linreg.h: catch situations where no solution is found
+
+2014-02-24  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin{1,2}D: Make minimal number of events per bin configurable
+
+2014-02-19  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin1D: Use weighted linear regression to determine D2
+
+2014-02-17  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin2D: Use weighted linear regression to determine D1
+
+2014-02-14  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin1D: Use weighted linear regression to determine D1
+
+2013-06-11  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Avoid extra copy of large input variables in C++ code
+
+2013-05-31  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries: be a little smarter to allow longer time series
+
+2013-03-26  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries: sampling frequency is a double now
+
+2013-03-15  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Switch documentation to in-source with roxygen2
+
+2012-12-17  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries: if no seed it given use a random one to preserve old behavior
+
+2012-12-10  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * timeseries: adding ability so set the seed for the RNG
+
+2012-11-20  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add David Bastine to package contributors
+
+2012-11-20  David Bastine <david.bastine@uni-oldenburg.de>
+    * timeseries: adding integration time step as optional argument
+
+2012-09-03  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin{1,2}D: adding option to set the number of threads
+
+2012-09-03  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * BugFix: Again: don't try to be smarter than you are in Langevin{1,2}D
+
+2012-08-24  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add function timeseries
+
+2012-08-15  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * BugFix: check for NAs when calculation mean_bins in Langevin2D
+
+2012-08-13  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * BugFix: Don't over-optimize in Langevin{1,2}D
+
+2012-08-09  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Langevin2D: Mayor rewrite of the function:
+        * Pure C++ now
+        * Output of D2 changed from (bins,bins,2,2) to (bins,bins,3)!
+
+2012-08-08  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * BugFix: corrected D2 estimation in Langevin1D
+    * Update documentation of Langevin{1,2}D
+
+2012-08-06  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Update documentation of Langevin1D
+
+2012-07-04  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Update the C++ and linker flags for Windows builds
+    * Byte-compile functions by default
+
+2012-05-25  Philip Rinn <philip.rinn@uni-oldenburg.de>
+    * Add generic functions for 1D and 2D Langevin analysis: Langevin{1,2}D
+    * Initiate the package

DESCRIPTION

+Package: Langevin
+Type: Package
+Title: Langevin Analysis in One and Two Dimensions
+Version: 1.1.1
+Date: 2015-11-03
+Authors@R: c(person('Philip', 'Rinn', email='philip.rinn@uni-oldenburg.de',
+    role=c('aut','cre')), person('Pedro G.', 'Lind', role='aut'),
+    person('David', 'Bastine', role='ctb'))
+Description: Estimate drift and diffusion functions from time series and
+    generate synthetic time series from given drift and diffusion coefficients.
+Encoding: UTF-8
+License: GPL (>= 2)
+LazyLoad: yes
+ByteCompile: yes
+NeedsCompilation: yes
+Depends:
+    R (>= 3.0.2)
+Imports:
+    Rcpp (>= 0.11.0)
+LinkingTo: Rcpp, RcppArmadillo (>= 0.4.600.0)
+RoxygenNote: 5.0.1

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+S3method(plot,Langevin)
+S3method(print,Langevin)
+S3method(summary,Langevin)
+export(Langevin1D)
+export(Langevin2D)
+export(timeseries1D)
+export(timeseries2D)
+import(Rcpp)
+importFrom(graphics,par)
+importFrom(graphics,plot)
+importFrom(stats,frequency)
+importFrom(stats,is.mts)
+importFrom(stats,is.ts)
+importFrom(stats,median)
+useDynLib(Langevin)

R/Langevin-package.r

+#' An \R package for stochastic data analysis
+#'
+#' The \pkg{Langevin} package provides functions to estimate drift and
+#' diffusion functions from data sets.
+#'
+#'
+#' This package was developed by the research group
+#' \emph{Turbulence, Wind energy and Stochastics} (TWiSt) at the Carl von
+#' Ossietzky University of Oldenburg (Germany).
+#'
+#' @section Mathematical Background: A wide range of dynamic systems can be
+#' described by a stochastic differential equation, the Langevin equation. The
+#' time derivative of the system trajectory \eqn{\dot{X}(t)} can be expressed as
+#' a sum of a deterministic part \eqn{D^{(1)}} and the product of a stochastic
+#' force \eqn{\Gamma(t)} and a weight coefficient \eqn{D^{(2)}}. The stochastic
+#' force \eqn{\Gamma(t)} is \eqn{\delta}-correlated Gaussian white noise.
+#'
+#' For stationary continuous Markov processes Siegert et al. and Friedrich et
+#' al. developed a method to reconstruct drift \eqn{D^{(1)}} and diffusion
+#' \eqn{D^{(2)}} directly from measured data.
+#'
+#' \deqn{ \dot{X}(t) = D^{(1)}(X(t),t) + \sqrt{D^{(2)}(X(t),t)}\,\Gamma(t)\quad
+#' \mathrm{with} } \deqn{ D^{(n)}(x,t) = \lim_{\tau \rightarrow 0}
+#' \frac{1}{\tau} M^{(n)}(x,t,\tau)\quad \mathrm{and} } \deqn{ M^{(n)}(x,t,\tau)
+#' = \frac{1}{n!} \langle (X(t+\tau) - x)^n \rangle |_{X(t) = x} }
+#'
+#' The Langevin equation should be interpreted in the way that for every time
+#' \eqn{t_i} where the system meets an arbitrary but fixed point \eqn{x} in
+#' phase space, \eqn{X(t_i+\tau)} is defined by the deterministic function
+#' \eqn{D^{(1)}(x)} and the stochastic function
+#' \eqn{\sqrt{D^{(2)}(x)}\Gamma(t_i)}. Both, \eqn{D^{(1)}(x)} and
+#' \eqn{D^{(2)}(x)} are constant for fixed \eqn{x}.
+#'
+#' One can integrate drift and diffusion numerically over small intervals. If
+#' the system is at time \eqn{t} in the state \eqn{x = X(t)} the drift can be
+#' calculated for small \eqn{\tau} by averaging over the difference of the
+#' system state at \eqn{t+\tau} and the state at \eqn{t}. The average has to be
+#' taken over the whole ensemble or in the stationary case over all \eqn{t =
+#' t_i} with \eqn{X(t_i) = x}. Diffusion can be calculated analogously.
+#'
+#'
+#' @name Langevin-package
+#' @aliases Langevin-package
+#' @docType package
+#' @author Philip Rinn
+#' @references
+#'
+#' \bold{A review of the Langevin method can be found at:}
+#'
+#' Friedrich, R.; et al. (2011) \emph{Approaching Complexity by Stochastic
+#' Methods: From Biological Systems to Turbulence}. Physics Reports, 506(5), 87–162.
+#'
+#' \bold{For further reading:}
+#'
+#' Risken, H. (1996) \emph{The Fokker-Planck equation}. Springer.
+#'
+#' Siegert, S.; et al. (1998) \emph{Analysis of data sets of stochastic
+#' systems}. Phys. Lett. A.
+#'
+#' Friedrich, R.; et al. (2000) \emph{Extracting model equations from
+#' experimental data}. Phys. Lett. A.
+NULL

R/Langevin1D.r

+#' Calculate the Drift and Diffusion of one-dimensional stochastic processes
+#'
+#' \code{Langevin1D} calculates the Drift and Diffusion vectors (with errors)
+#' for a given time series.
+#'
+#'
+#' @param data a vector containing the time series or a time-series object.
+#' @param bins a scalar denoting the number of \code{bins} to calculate the
+#' conditional moments on.
+#' @param steps a vector giving the \eqn{\tau} steps to calculate the
+#' conditional moments (in samples (=\eqn{\tau * sf})).
+#' @param sf a scalar denoting the sampling frequency (optional if \code{data}
+#' is a time-series object).
+#' @param bin_min a scalar denoting the minimal number of events per \code{bin}.
+#' Defaults to \code{100}.
+#' @param reqThreads a scalar denoting how many threads to use. Defaults to
+#' \code{-1} which means all available cores.
+#'
+#' @return \code{Langevin1D} returns a list with thirteen components:
+#' @return \item{D1}{a vector of the Drift coefficient for each \code{bin}.}
+#' @return \item{eD1}{a vector of the error of the Drift coefficient for each
+#' \code{bin}.}
+#' @return \item{D2}{a vector of the Diffusion coefficient for each \code{bin}.}
+#' @return \item{eD2}{a vector of the error of the Driffusion coefficient for
+#' each \code{bin}.}
+#' @return \item{D4}{a vector of the fourth Kramers-Moyal coefficient for each
+#' \code{bin}.}
+#' @return \item{mean_bin}{a vector of the mean value per \code{bin}.}
+#' @return \item{density}{a vector of the number of events per \code{bin}.}
+#' @return \item{M1}{a matrix of the first conditional moment for each
+#' \eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
+#' @return \item{eM1}{a matrix of the error of the first conditional moment
+#' for each \eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
+#' @return \item{M2}{a matrix of the second conditional moment for each
+#' \eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
+#' @return \item{eM2}{a matrix of the error of the second conditional moment
+#' for each \eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
+#' @return \item{M4}{a matrix of the fourth conditional moment for each
+#' \eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
+#' @return \item{U}{a vector of the \code{bin} borders.}
+#'
+#' @author Philip Rinn
+#' @seealso \code{\link{Langevin2D}}
+#' @examples
+#'
+#' # Set number of bins, steps and the sampling frequency
+#' bins <- 20;
+#' steps <- c(1:5);
+#' sf <- 1000;
+#'
+#' #### Linear drift, constant diffusion
+#'
+#' # Generate a time series with linear D^1 = -x and constant D^2 = 1
+#' x <- timeseries1D(N=1e6, d11=-1, d20=1, sf=sf);
+#' # Do the analysis
+#' est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
+#' # Plot the result and add the theoretical expectation as red line
+#' plot(est$mean_bin, est$D1);
+#' lines(est$mean_bin, -est$mean_bin, col='red');
+#' plot(est$mean_bin, est$D2);
+#' abline(h=1, col='red');
+#'
+#' #### Cubic drift, constant diffusion
+#'
+#' # Generate a time series with cubic D^1 = x - x^3 and constant D^2 = 1
+#' x <- timeseries1D(N=1e6, d13=-1, d11=1, d20=1, sf=sf);
+#' # Do the analysis
+#' est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
+#' # Plot the result and add the theoretical expectation as red line
+#' plot(est$mean_bin, est$D1);
+#' lines(est$mean_bin, est$mean_bin - est$mean_bin^3, col='red');
+#' plot(est$mean_bin, est$D2);
+#' abline(h=1, col='red');
+#' @import Rcpp
+#' @importFrom stats frequency is.ts
+#' @useDynLib Langevin
+#' @export
+Langevin1D <- function(data, bins, steps,
+                       sf=ifelse(is.ts(data), frequency(data), 1), bin_min=100,
+                       reqThreads=-1) {
+    .Call('Langevin_Langevin1D', PACKAGE='Langevin', data, bins, steps, sf,
+          bin_min, reqThreads)
+}

R/Langevin2D.r

+#' Calculate the Drift and Diffusion of two-dimensional stochastic processes
+#'
+#' \code{Langevin2D} calculates the Drift (with error) and Diffusion matrices
+#' for given time series.
+#'
+#'
+#' @param data a matrix containing the time series as columns or a time-series
+#' object.
+#' @param bins a scalar denoting the number of \code{bins} to calculate Drift
+#' and Diffusion on.
+#' @param steps a vector giving the \eqn{\tau} steps to calculate the moments
+#' (in samples).
+#' @param sf a scalar denoting the sampling frequency (optional if \code{data}
+#' is a time-series object).
+#' @param bin_min a scalar denoting the minimal number of events per \code{bin}.
+#' Defaults to \code{100}.
+#' @param reqThreads a scalar denoting how many threads to use. Defaults to
+#' \code{-1} which means all available cores.
+#'
+#' @return \code{Langevin2D} returns a list with nine components:
+#' @return \item{D1}{a tensor with all values of the drift coefficient.
+#' Dimension is \code{bins} x \code{bins} x 2. The first
+#' \code{bins} x \code{bins} elements define the drift \eqn{D^{(1)}_{1}}
+#' for the first variable and the rest define the drift \eqn{D^{(1)}_{2}}
+#' for the second variable.}
+#' @return \item{eD1}{a tensor with all estimated errors of the drift
+#' coefficient. Dimension is \code{bins} x \code{bins} x 2. Same expression as
+#' above.}
+#' @return \item{D2}{a tensor with all values of the diffusion coefficient.
+#' Dimension is \code{bins} x \code{bins} x 3. The first
+#' \code{bins} x \code{bins} elements define the diffusion \eqn{D^{(2)}_{11}},
+#' the second \code{bins} x \code{bins} elements define the diffusion
+#' \eqn{D^{(2)}_{22}} and the rest define the diffusion
+#' \eqn{D^{(2)}_{12} = D^{(2)}_{21}}.}
+#' @return \item{eD2}{a tensor with all estimated errors of the driffusion
+#' coefficient. Dimension is \code{bins} x \code{bins} x 3. Same expression as
+#' above.}
+#' @return \item{mean_bin}{a matrix of the mean value per \code{bin}.
+#' Dimension is \code{bins} x \code{bins} x 2. The first
+#' \code{bins} x \code{bins} elements define the mean for the first variable
+#' and the rest for the second variable.}
+#' @return \item{density}{a matrix of the number of events per \code{bin}.
+#' Rows label the \code{bin} of the first variable and columns the second
+#' variable.}
+#' @return \item{M1}{a tensor of the first moment for each \code{bin} (line
+#' label) and  each \eqn{\tau} step (column label). Dimension is
+#' \code{bins} x \code{bins} x 2\code{length(steps)}.}
+#' @return \item{eM1}{a tensor of the standard deviation of the first
+#' moment for each bin (line label) and  each \eqn{\tau} step (column label).
+#' Dimension is \code{bins} x \code{bins} x 2\code{length(steps)}.}
+#' @return \item{M2}{a tensor of the second moment for each bin (line
+#' label) and  each \eqn{\tau} step (column label). Dimension is
+#' \code{bins} x \code{bins} x 3\code{length(steps)}.}
+#' @return \item{U}{a matrix of the \code{bin} borders}
+#'
+#' @author Philip Rinn
+#' @seealso \code{\link{Langevin1D}}
+#' @import Rcpp
+#' @importFrom stats frequency is.mts
+#' @useDynLib Langevin
+#' @export
+Langevin2D <- function(data, bins, steps,
+                       sf=ifelse(is.mts(data), frequency(data), 1), bin_min=100,
+                       reqThreads=-1) {
+    if(nrow(data) == 2)
+        stop("Time series have to be arranged in colums now. Please adopt your code!")
+    .Call('Langevin_Langevin2D', PACKAGE='Langevin', data, bins, steps, sf,
+          bin_min, reqThreads)
+}

R/RcppExports.R

+# This file was generated by Rcpp::compileAttributes
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+.Langevin1D <- function(data, bins, steps, sf, bin_min, reqThreads) {
+    .Call('Langevin_Langevin1D', PACKAGE = 'Langevin', data, bins, steps, sf, bin_min, reqThreads)
+}
+
+.Langevin2D <- function(data, bins, steps, sf, bin_min, reqThreads) {
+    .Call('Langevin_Langevin2D', PACKAGE = 'Langevin', data, bins, steps, sf, bin_min, reqThreads)
+}
+
+.timeseries1D <- function(N, startpoint, d13, d12, d11, d10, d22, d21, d20, sf, dt) {
+    .Call('Langevin_timeseries1D', PACKAGE = 'Langevin', N, startpoint, d13, d12, d11, d10, d22, d21, d20, sf, dt)
+}
+
+.timeseries2D <- function(N, startpointx, startpointy, D1_1, D1_2, g_11, g_12, g_21, g_22, sf, dt) {
+    .Call('Langevin_timeseries2D', PACKAGE = 'Langevin', N, startpointx, startpointy, D1_1, D1_2, g_11, g_12, g_21, g_22, sf, dt)
+}
+

R/plot.Langevin.r

+#' Plot estimated drift and diffusion coefficients
+#'
+#' plot method for class "Langevin".
+#' This method is only implemented for one-dimensional analysis for now.
+#'
+#'
+#' @param x an object of class "Langevin".
+#' @param pch Either an integer specifying a symbol or a single character to be
+#' used as the default in plotting points. See \link{points} for possible values
+#' and their interpretation. Default is \code{pch = 20}.
+#' @param ... Arguments to be passed to methods, such as \code{\link{par}}.
+#'
+#'
+#' @author Philip Rinn
+#' @importFrom graphics par plot
+#' @export
+plot.Langevin <- function(x, pch=20, ...) {
+    if(dim(x$D1)[2] > 1)
+        stop("Plotting is only implemented for the one-dimensional Langevin
+             Approach")
+    orig_mar = par("mar")
+    orig_mfrow = par("mfrow")
+    par(mar=c(4.2, 5.4, 0.5, 0.5))
+    par(mfrow=c(1,2))
+    plot(x$mean_bin, x$D1, xlab="x [a.u.]",
+         ylab=expression(paste("Drift coefficient ",D^(1), "(x) [",1/s,"]")),
+         pch=pch, ...)
+    plot(x$mean_bin, x$D2, xlab="x [a.u.]",
+           ylab=expression(paste("Diffusion coefficient ",D^(2),
+                                 "(x) [",1/s^2,"]")), pch=pch, ...)
+    par(mar=orig_mar)
+    par(mfrow=orig_mfrow)
+}

R/print.Langevin.r

+#' Print estimated drift and diffusion coefficients
+#'
+#' print method for class "Langevin".
+#'
+#'
+#' @param x an object of class "Langevin".
+#' @param digits integer, used for number formatting with \code{\link{signif}()}.
+#' @param ... further arguments to be passed to or from other methods. They are
+#' ignored in this function.
+#'
+#' @return The function \code{print.Langevin()} returns an overview of the
+#' estimated drift and diffusion coefficients.
+#'
+#' @author Philip Rinn
+#' @export
+print.Langevin <- function(x, digits=max(3, getOption("digits")-3), ...) {
+    cat(paste0("Drift and diffusion estimates for the ",
+               ifelse(dim(x$D1)[2] > 1, "two", "one"),
+               " dimensional Langevin Approach\n"),
+        paste0("Number of bins: ", dim(x$D1)[1],
+               if(dim(x$D1)[2] > 1) paste0("x", dim(x$D1)[2]), "\n"),
+        paste0("Range of the bins: ", signif(range(x$U,na.rm=T)[1], digits),
+               " ... ", signif(range(x$U,na.rm=T)[2], digits),"\n"),
+        paste0("Range of D1: ", signif(range(x$D1,na.rm=T)[1], digits),
+               " ... ", signif(range(x$D1,na.rm=T)[2], digits),"\n"),
+        paste0("Range of D2: ", signif(range(x$D2,na.rm=T)[1], digits),
+               " ... ", signif(range(x$D2,na.rm=T)[2], digits),"\n")
+        )
+}

R/summary.Langevin.r

+#' Summarize estimated drift and diffusion coefficients
+#'
+#' summary method for class "Langevin".
+#'
+#'
+#' @param object an object of class "Langevin".
+#' @param ... arguments to be passed to or from other methods. They are
+#' ignored in this function.
+#' @param digits integer, used for number formatting with \code{\link{signif}()}.
+#'
+#' @return The function \code{summary.Langevin()} returns a sumamry of the
+#' estimated drift and diffusion coefficients
+#'
+#' @author Philip Rinn
+#' @importFrom stats median
+#' @export
+summary.Langevin <- function(object, ..., digits=max(3, getOption("digits") - 3)) {
+    cat(paste0(" Number of bins: ", dim(object$D1)[1],
+               if(dim(object$D1)[2] > 1) paste0("x", dim(object$D1)[2]), "\n"),
+        paste0("Population of the bins:\n",
+              "\tMin.  : ", min(object$density, na.rm=T), "\n",
+              "\tMedian: ", round(median(object$density, na.rm=T)), "\n",
+              "\tMean  : ", round(mean(object$density, na.rm=T)), "\n",
+              "\tMax.  : ", max(object$density, na.rm=T), "\n"),
+        paste0("Number of NA's for D1: ", length(which(is.na(object$D1))),
+               "\n"),
+        paste0("Number of NA's for D2: ", length(which(is.na(object$D2))),
+               "\n"),
+        if(dim(object$D1)[2] == 1) paste0("Ratio between D4 and D2^2:\n",
+              "\tMin.  : ", signif(min(object$D4/object$D2^2, na.rm=T), digits), "\n",
+              "\tMedian: ", signif(median(object$D4/object$D2^2, na.rm=T), digits), "\n",
+              "\tMean  : ", signif(mean(object$D4/object$D2^2, na.rm=T), digits), "\n",
+              "\tMax.  : ", signif(max(object$D4/object$D2^2, na.rm=T), digits), "\n")
+        )
+}

R/timeseries1D.r

+#' Generate a 1D Langevin process
+#'
+#' \code{timeseries1D} generates a one-dimensional Langevin process using a
+#' simple Euler integration. The drift function is a cubic polynomial, the
+#' diffusion funcation a quadratic.
+#'
+#'
+#' @param N a scalar denoting the length of the time-series to generate.
+#' @param startpoint a scalar denoting the starting point of the time series.
+#' @param d13,d12,d11,d10 scalars denoting the coefficients for the drift polynomial.
+#' @param d22,d21,d20 scalars denoting the coefficients for the diffusion polynomial.
+#' @param sf a scalar denoting the sampling frequency.
+#' @param dt a scalar denoting the maximal time step of integration. Default
+#' \code{dt=0} yields \code{dt=1/sf}.
+#'
+#' @return \code{timeseries1D} returns a time-series object of length
+#' \code{N} with the generated time-series.
+#'
+#' @author Philip Rinn
+#' @seealso \code{\link{timeseries2D}}
+#' @examples
+#' # Generate standardized Ornstein-Uhlenbeck-Process (d11=-1, d20=1)
+#' # with integration time step 0.01 and sampling frequency 1
+#' s <- timeseries1D(N=1e4, sf=1, dt=0.01);
+#' t <- 1:1e4;
+#' plot(t, s, t="l", main=paste("mean:", mean(s), " var:", var(s)));
+#' @import Rcpp
+#' @useDynLib Langevin
+#' @export
+timeseries1D <- function(N, startpoint = 0, d13 = 0, d12 = 0, d11 = -1, d10 = 0,
+                         d22 = 0, d21 = 0, d20 = 1, sf = 1000, dt = 0) {
+    .Call('Langevin_timeseries1D', PACKAGE = 'Langevin', N, startpoint, d13,
+          d12, d11, d10, d22, d21, d20, sf, dt)
+}

R/timeseries2D.r

+# Copyright (C) 2015  Philip Rinn
+# Copyright (C) 2015  Carl von Ossietzy Universität Oldenburg
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, see <https://www.gnu.org/licenses/gpl-2.0>.
+
+#' Generate a 2D Langevin process
+#'
+#' \code{timeseries2D} generates a two-dimensional Langevin process using a
+#' simple Euler integration. The drift function is a cubic polynomial, the
+#' diffusion function a quadratic.
+#'
+#' The elements \eqn{a_{ij}} of the matrices are defined by the corresponding
+#' equations for the drift and diffusion terms:
+#'
+#' \deqn{D^1_{1,2} = \sum_{i,j=1}^4 a_{ij} x_1^{(i-1)}x_2^{(j-1)} }
+#'
+#' with \eqn{a_{ij} = 0} for \eqn{ i + j > 5}.
+#'
+#' \deqn{g_{11,12,21,22} = \sum_{i,j=1}^3 a_{ij} x_1^{(i-1)}x_2^{(j-1)} }
+#'
+#' with \eqn{a_{ij} = 0} for \eqn{ i + j > 4}
+#'
+#' @param N a scalar denoting the length of the time-series to generate.
+#' @param startpointx a scalar denoting the starting point of the time series x.
+#' @param startpointy a scalar denoting the starting point of the time series y.
+#' @param D1_1 a 4x4 matrix denoting the coefficients of D1 for x.
+#' @param D1_2 a 4x4 matrix denoting the coefficients of D1 for y.
+#' @param g_11 a 3x3 matrix denoting the coefficients of g11 for x.
+#' @param g_12 a 3x3 matrix denoting the coefficients of g12 for x.
+#' @param g_21 a 3x3 matrix denoting the coefficients of g21 for y.
+#' @param g_22 a 3x3 matrix denoting the coefficients of g22 for y.
+#' @param sf a scalar denoting the sampling frequency.
+#' @param dt a scalar denoting the maximal time step of integration. Default
+#' \code{dt=0} yields \code{dt=1/sf}.
+#'
+#' @return \code{timeseries2D} returns a time-series object with the generated
+#' time-series as colums.
+#'
+#' @author Philip Rinn
+#' @seealso \code{\link{timeseries1D}}
+#' @import Rcpp
+#' @useDynLib Langevin
+#' @export
+timeseries2D <- function(N, startpointx=0, startpointy=0,
+                         D1_1=matrix(c(0,-1,rep(0,14)),nrow=4),
+                         D1_2=matrix(c(0,0,0,0,-1,rep(0,11)),nrow=4),
+                         g_11=matrix(c(1,0,0,0,0,0,0,0,0),nrow=3),
+                         g_12=matrix(c(0,0,0,0,0,0,0,0,0),nrow=3),
+                         g_21=matrix(c(0,0,0,0,0,0,0,0,0),nrow=3),
+                         g_22=matrix(c(1,0,0,0,0,0,0,0,0),nrow=3),
+                         sf=1000, dt=0) {
+    .Call('Langevin_timeseries2D', PACKAGE = 'Langevin', N, startpointx,
+          startpointy, D1_1, D1_2, g_11, g_12, g_21, g_22, sf, dt)
+}

cleanup

+#!/bin/sh
+
+rm -f  src/*.o src/*.so
+

inst/CITATION

+year <- sub("-.*", "", meta$Date)
+note <- sprintf("R package version %s", meta$Version)
+
+bibentry(
+    header = "To cite the 'Langevin' package in publications use:",
+    bibtype = "Manual",
+    key = "Rinn2015",
+    title = "Langevin Analysis in One and Two Dimensions",
+    author  = c(person('Philip', 'Rinn', role=c('aut','cre')),
+                person('Pedro G.', 'Lind', role='aut'),
+                person('David', 'Bastine', role='ctb')),
+    year = year,
+    note = note,
+    url = "https://CRAN.R-project.org/package=Langevin")
+
+bibentry(
+    header = "To cite the 'Langevin Approach' in publications use:",
+    bibtype = "Article",
+    key = "Friedrich2011",
+    title = "Approaching complexity by stochastic methods: From biological systems to turbulence",
+    author  = c(person("Rudolf", "Friedrich"),
+                person("Joachim", "Peinke"),
+                person("Muhammad", "Sahimi"),
+                person("Mohammed", "Reza Rahimi Tabar")),
+    journal = "Physics Reports",
+    year = 2011,
+    number = 5,
+    pages = "87--162",
+    volume = 506,
+    doi = "10.1016/j.physrep.2011.05.003",
+    issn = "0370-1573",
+    publisher = "Elsevier BV")
+

man/Langevin-package.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Langevin-package.r
+\docType{package}
+\name{Langevin-package}
+\alias{Langevin-package}
+\title{An \R package for stochastic data analysis}
+\description{
+The \pkg{Langevin} package provides functions to estimate drift and
+diffusion functions from data sets.
+}
+\details{
+This package was developed by the research group
+\emph{Turbulence, Wind energy and Stochastics} (TWiSt) at the Carl von
+Ossietzky University of Oldenburg (Germany).
+}
+\section{Mathematical Background}{
+ A wide range of dynamic systems can be
+described by a stochastic differential equation, the Langevin equation. The
+time derivative of the system trajectory \eqn{\dot{X}(t)} can be expressed as
+a sum of a deterministic part \eqn{D^{(1)}} and the product of a stochastic
+force \eqn{\Gamma(t)} and a weight coefficient \eqn{D^{(2)}}. The stochastic
+force \eqn{\Gamma(t)} is \eqn{\delta}-correlated Gaussian white noise.
+
+For stationary continuous Markov processes Siegert et al. and Friedrich et
+al. developed a method to reconstruct drift \eqn{D^{(1)}} and diffusion
+\eqn{D^{(2)}} directly from measured data.
+
+\deqn{ \dot{X}(t) = D^{(1)}(X(t),t) + \sqrt{D^{(2)}(X(t),t)}\,\Gamma(t)\quad
+\mathrm{with} } \deqn{ D^{(n)}(x,t) = \lim_{\tau \rightarrow 0}
+\frac{1}{\tau} M^{(n)}(x,t,\tau)\quad \mathrm{and} } \deqn{ M^{(n)}(x,t,\tau)
+= \frac{1}{n!} \langle (X(t+\tau) - x)^n \rangle |_{X(t) = x} }
+
+The Langevin equation should be interpreted in the way that for every time
+\eqn{t_i} where the system meets an arbitrary but fixed point \eqn{x} in
+phase space, \eqn{X(t_i+\tau)} is defined by the deterministic function
+\eqn{D^{(1)}(x)} and the stochastic function
+\eqn{\sqrt{D^{(2)}(x)}\Gamma(t_i)}. Both, \eqn{D^{(1)}(x)} and
+\eqn{D^{(2)}(x)} are constant for fixed \eqn{x}.
+
+One can integrate drift and diffusion numerically over small intervals. If
+the system is at time \eqn{t} in the state \eqn{x = X(t)} the drift can be
+calculated for small \eqn{\tau} by averaging over the difference of the
+system state at \eqn{t+\tau} and the state at \eqn{t}. The average has to be
+taken over the whole ensemble or in the stationary case over all \eqn{t =
+t_i} with \eqn{X(t_i) = x}. Diffusion can be calculated analogously.