Commit f8751b9a authored by Philip Rinn's avatar Philip Rinn
Browse files

Langevin1D: add option to use the kernel based Nadaraya-Watson estimator

parent dfc78962
2021-10-18 Philip Rinn <philip.rinn@uni-oldenburg.de>
* Tag as version 1.3.0
* Langevin1D: add option to use the kernel based Nadaraya-Watson estimator
2021-09-04 Philip Rinn <philip.rinn@uni-oldenburg.de>
* Fix documentation: eD2 is not implemented for Langevin2D
* Fix spelling errors
2019-01-02 Philip Rinn <philip.rinn@uni-oldenburg.de>
* Tag as version 1.2.1
* Remove calls to `Rcpp:::CxxFlags()` and `Rcpp:::LdFlags()`
......
Package: Langevin
Type: Package
Title: Langevin Analysis in One and Two Dimensions
Version: 1.2.1
Date: 2019-01-02
Version: 1.3.0
Date: 2021-10-18
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'))
......@@ -11,7 +11,6 @@ Description: Estimate drift and diffusion functions from time series and
Encoding: UTF-8
License: GPL (>= 2)
URL: https://gitlab.uni-oldenburg.de/TWiSt/Langevin
BugReports: https://gitlab.uni-oldenburg.de/TWiSt/Langevin/issues
LazyLoad: yes
ByteCompile: yes
NeedsCompilation: yes
......@@ -20,4 +19,4 @@ Depends:
Imports:
Rcpp (>= 0.11.0)
LinkingTo: Rcpp, RcppArmadillo (>= 0.4.600.0)
RoxygenNote: 6.1.1
RoxygenNote: 7.1.2
......@@ -10,6 +10,7 @@ export(timeseries2D)
import(Rcpp)
importFrom(graphics,par)
importFrom(graphics,plot)
importFrom(stats,bw.nrd)
importFrom(stats,frequency)
importFrom(stats,is.mts)
importFrom(stats,is.ts)
......
......@@ -59,4 +59,7 @@
#'
#' Friedrich, R.; et al. (2000) \emph{Extracting model equations from
#' experimental data}. Phys. Lett. A.
#'
#' Honisch, C.; Friedrich, R. (2011). \emph{Estimation of Kramers-Moyal
#' coefficients at low sampling rates.}. Physical Review E, 83(6), 066701.
NULL
......@@ -8,15 +8,23 @@
#' @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})).
#' conditional moments (in samples (=\eqn{\tau * sf})). Only used if
#' \code{kernel} is \code{FALSE}.
#' @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.
#' \code{-1} which means all available cores. Only used if \code{kernel} is
#' \code{FALSE}.
#' @param kernel a logical denoting if the kernel based Nadaraya-Watson
#' estimator should be used to calculate drift and diffusion vectors.
#' @param h a scalar denoting the bandwidth of the data. Defaults to Scott's
#' variation of Silverman's rule of thumb. Only used if \code{kernel} is
#' \code{TRUE}.
#'
#' @return \code{Langevin1D} returns a list with thirteen components:
#' @return \code{Langevin1D} returns a list with thirteen (six if \code{kernel}
#' is \code{TRUE}) 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}.}
......@@ -26,57 +34,73 @@
#' @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{density}{a vector of the number of events per \code{bin}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{M1}{a matrix of the first conditional moment for each
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.}
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{eM1}{a matrix of the error of the first conditional moment
#' for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.}
#' for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{M2}{a matrix of the second conditional moment for each
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.}
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{eM2}{a matrix of the error of the second conditional moment
#' for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.}
#' for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{M4}{a matrix of the fourth conditional moment for each
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.}
#' @return \item{U}{a vector of the \code{bin} borders.}
#' \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
#' If \code{kernel} is \code{FALSE}.}
#' @return \item{U}{a vector of the \code{bin} borders.
#' If \code{kernel} is \code{FALSE}.}
#'
#' @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;
#' 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);
#' x <- timeseries1D(N = 1e6, d11 = -1, d20 = 1, sf = sf)
#' # Do the analysis
#' est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
#' est <- Langevin1D(data = x, bins = bins, steps = steps, sf = sf)
#' # 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');
#' 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);
#' x <- timeseries1D(N = 1e6, d13 = -1, d11 = 1, d20 = 1, sf = sf)
#' # Do the analysis
#' est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
#' est <- Langevin1D(data = x, bins = bins, steps = steps, sf = sf)
#' # 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');
#' 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
#' @importFrom stats frequency is.ts bw.nrd
#' @useDynLib Langevin, .registration=TRUE
#' @export
Langevin1D <- function(data, bins, steps,
sf=ifelse(is.ts(data), frequency(data), 1), bin_min=100,
reqThreads=-1) {
.Call('_Langevin_Langevin1D', data, bins, steps, sf, bin_min, reqThreads)
reqThreads=-1, kernel=FALSE, h) {
if (kernel) {
steps <- NULL
if (missing(h)) {
h <- bw.nrd(data)
}
.Call("_Langevin_kernel1D", data, bins, sf, h)
} else {
h <- NULL
.Call("_Langevin_Langevin1D", data, bins, steps, sf, bin_min, reqThreads)
}
}
......@@ -9,6 +9,10 @@
.Call(`_Langevin_Langevin2D`, data, bins, steps, sf, bin_min, reqThreads)
}
.kernel1D <- function(data, bins, sf, h) {
.Call(`_Langevin_kernel1D`, data, bins, sf, h)
}
.timeseries1D <- function(N, startpoint, d13, d12, d11, d10, d22, d21, d20, sf, dt) {
.Call(`_Langevin_timeseries1D`, N, startpoint, d13, d12, d11, d10, d22, d21, d20, sf, dt)
}
......
......@@ -19,8 +19,8 @@ print.Langevin <- function(x, digits=max(3, getOption("digits")-3), ...) {
" 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 the bins: ", ifelse("U" %in% names(x), signif(range(x$U,na.rm=T)[1], digits), signif(range(x$mean_bin,na.rm=T)[1], digits)),
" ... ", ifelse("U" %in% names(x), signif(range(x$U,na.rm=T)[2], digits), signif(range(x$mean_bin,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),
......
......@@ -17,11 +17,12 @@
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"),
if("density" %in% names(object)) {
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"),
"\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))),
......
......@@ -60,6 +60,9 @@ systems}. Phys. Lett. A.
Friedrich, R.; et al. (2000) \emph{Extracting model equations from
experimental data}. Phys. Lett. A.
Honisch, C.; Friedrich, R. (2011). \emph{Estimation of Kramers-Moyal
coefficients at low sampling rates.}. Physical Review E, 83(6), 066701.
}
\author{
Philip Rinn
......
......@@ -4,8 +4,16 @@
\alias{Langevin1D}
\title{Calculate the Drift and Diffusion of one-dimensional stochastic processes}
\usage{
Langevin1D(data, bins, steps, sf = ifelse(is.ts(data), frequency(data),
1), bin_min = 100, reqThreads = -1)
Langevin1D(
data,
bins,
steps,
sf = ifelse(is.ts(data), frequency(data), 1),
bin_min = 100,
reqThreads = -1,
kernel = FALSE,
h
)
}
\arguments{
\item{data}{a vector containing the time series or a time-series object.}
......@@ -14,7 +22,8 @@ Langevin1D(data, bins, steps, sf = ifelse(is.ts(data), frequency(data),
conditional moments on.}
\item{steps}{a vector giving the \eqn{\tau} steps to calculate the
conditional moments (in samples (=\eqn{\tau * sf})).}
conditional moments (in samples (=\eqn{\tau * sf})). Only used if
\code{kernel} is \code{FALSE}.}
\item{sf}{a scalar denoting the sampling frequency (optional if \code{data}
is a time-series object).}
......@@ -23,10 +32,19 @@ is a time-series object).}
Defaults to \code{100}.}
\item{reqThreads}{a scalar denoting how many threads to use. Defaults to
\code{-1} which means all available cores.}
\code{-1} which means all available cores. Only used if \code{kernel} is
\code{FALSE}.}
\item{kernel}{a logical denoting if the kernel based Nadaraya-Watson
estimator should be used to calculate drift and diffusion vectors.}
\item{h}{a scalar denoting the bandwidth of the data. Defaults to Scott's
variation of Silverman's rule of thumb. Only used if \code{kernel} is
\code{TRUE}.}
}
\value{
\code{Langevin1D} returns a list with thirteen components:
\code{Langevin1D} returns a list with thirteen (six if \code{kernel}
is \code{TRUE}) components:
\item{D1}{a vector of the Drift coefficient for each \code{bin}.}
......@@ -35,7 +53,7 @@ Defaults to \code{100}.}
\item{D2}{a vector of the Diffusion coefficient for each \code{bin}.}
\item{eD2}{a vector of the error of the Driffusion coefficient for
\item{eD2}{a vector of the error of the Diffusion coefficient for
each \code{bin}.}
\item{D4}{a vector of the fourth Kramers-Moyal coefficient for each
......@@ -43,24 +61,31 @@ each \code{bin}.}
\item{mean_bin}{a vector of the mean value per \code{bin}.}
\item{density}{a vector of the number of events per \code{bin}.}
\item{density}{a vector of the number of events per \code{bin}.
If \code{kernel} is \code{FALSE}.}
\item{M1}{a matrix of the first conditional moment for each
\eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
\eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
If \code{kernel} is \code{FALSE}.}
\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}.}
for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
If \code{kernel} is \code{FALSE}.}
\item{M2}{a matrix of the second conditional moment for each
\eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
\eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
If \code{kernel} is \code{FALSE}.}
\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}.}
for each \eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
If \code{kernel} is \code{FALSE}.}
\item{M4}{a matrix of the fourth conditional moment for each
\eqn{\tau}. Rows corespond to \code{bin}, columns to \eqn{\tau}.}
\eqn{\tau}. Rows correspond to \code{bin}, columns to \eqn{\tau}.
If \code{kernel} is \code{FALSE}.}
\item{U}{a vector of the \code{bin} borders.}
\item{U}{a vector of the \code{bin} borders.
If \code{kernel} is \code{FALSE}.}
}
\description{
\code{Langevin1D} calculates the Drift and Diffusion vectors (with errors)
......@@ -69,33 +94,33 @@ for a given time series.
\examples{
# Set number of bins, steps and the sampling frequency
bins <- 20;
steps <- c(1:5);
sf <- 1000;
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);
x <- timeseries1D(N = 1e6, d11 = -1, d20 = 1, sf = sf)
# Do the analysis
est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
est <- Langevin1D(data = x, bins = bins, steps = steps, sf = sf)
# 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');
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);
x <- timeseries1D(N = 1e6, d13 = -1, d11 = 1, d20 = 1, sf = sf)
# Do the analysis
est <- Langevin1D(x, bins, steps, sf, reqThreads=2);
est <- Langevin1D(data = x, bins = bins, steps = steps, sf = sf)
# 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');
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")
}
\seealso{
\code{\link{Langevin2D}}
......
......@@ -4,8 +4,14 @@
\alias{Langevin2D}
\title{Calculate the Drift and Diffusion of two-dimensional stochastic processes}
\usage{
Langevin2D(data, bins, steps, sf = ifelse(is.mts(data), frequency(data),
1), bin_min = 100, reqThreads = -1)
Langevin2D(
data,
bins,
steps,
sf = ifelse(is.mts(data), frequency(data), 1),
bin_min = 100,
reqThreads = -1
)
}
\arguments{
\item{data}{a matrix containing the time series as columns or a time-series
......@@ -46,10 +52,6 @@ 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}}.}
\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.}
\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
......
......@@ -4,8 +4,7 @@
\alias{print.Langevin}
\title{Print estimated drift and diffusion coefficients}
\usage{
\method{print}{Langevin}(x, digits = max(3, getOption("digits") - 3),
...)
\method{print}{Langevin}(x, digits = max(3, getOption("digits") - 3), ...)
}
\arguments{
\item{x}{an object of class "Langevin".}
......
......@@ -4,8 +4,7 @@
\alias{summary.Langevin}
\title{Summarize estimated drift and diffusion coefficients}
\usage{
\method{summary}{Langevin}(object, ..., digits = max(3,
getOption("digits") - 3))
\method{summary}{Langevin}(object, ..., digits = max(3, getOption("digits") - 3))
}
\arguments{
\item{object}{an object of class "Langevin".}
......@@ -16,7 +15,7 @@ ignored in this function.}
\item{digits}{integer, used for number formatting with \code{\link{signif}()}.}
}
\value{
The function \code{summary.Langevin()} returns a sumamry of the
The function \code{summary.Langevin()} returns a summary of the
estimated drift and diffusion coefficients
}
\description{
......
......@@ -4,8 +4,19 @@
\alias{timeseries1D}
\title{Generate a 1D Langevin process}
\usage{
timeseries1D(N, startpoint = 0, d13 = 0, d12 = 0, d11 = -1,
d10 = 0, d22 = 0, d21 = 0, d20 = 1, sf = 1000, dt = 0)
timeseries1D(
N,
startpoint = 0,
d13 = 0,
d12 = 0,
d11 = -1,
d10 = 0,
d22 = 0,
d21 = 0,
d20 = 1,
sf = 1000,
dt = 0
)
}
\arguments{
\item{N}{a scalar denoting the length of the time-series to generate.}
......@@ -28,7 +39,7 @@ timeseries1D(N, startpoint = 0, d13 = 0, d12 = 0, d11 = -1,
\description{
\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.
diffusion function a quadratic.
}
\examples{
# Generate standardized Ornstein-Uhlenbeck-Process (d11=-1, d20=1)
......
......@@ -4,13 +4,19 @@
\alias{timeseries2D}
\title{Generate a 2D Langevin process}
\usage{
timeseries2D(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),
timeseries2D(
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)
g_22 = matrix(c(1, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 3),
sf = 1000,
dt = 0
)
}
\arguments{
\item{N}{a scalar denoting the length of the time-series to generate.}
......@@ -38,7 +44,7 @@ timeseries2D(N, startpointx = 0, startpointy = 0, D1_1 = matrix(c(0,
}
\value{
\code{timeseries2D} returns a time-series object with the generated
time-series as colums.
time-series as columns.
}
\description{
\code{timeseries2D} generates a two-dimensional Langevin process using a
......
......@@ -6,6 +6,11 @@
using namespace Rcpp;
#ifdef RCPP_USE_GLOBAL_ROSTREAM
Rcpp::Rostream<true>& Rcpp::Rcout = Rcpp::Rcpp_cout_get();
Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
#endif
// Langevin1D
List Langevin1D(const arma::vec& data, const int& bins, const arma::vec& steps, const double& sf, const int& bin_min, int reqThreads);
RcppExport SEXP _Langevin_Langevin1D(SEXP dataSEXP, SEXP binsSEXP, SEXP stepsSEXP, SEXP sfSEXP, SEXP bin_minSEXP, SEXP reqThreadsSEXP) {
......@@ -38,6 +43,20 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
// kernel1D
List kernel1D(const arma::vec& data, const int& bins, const double& sf, const double& h);
RcppExport SEXP _Langevin_kernel1D(SEXP dataSEXP, SEXP binsSEXP, SEXP sfSEXP, SEXP hSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const arma::vec& >::type data(dataSEXP);
Rcpp::traits::input_parameter< const int& >::type bins(binsSEXP);
Rcpp::traits::input_parameter< const double& >::type sf(sfSEXP);
Rcpp::traits::input_parameter< const double& >::type h(hSEXP);
rcpp_result_gen = Rcpp::wrap(kernel1D(data, bins, sf, h));
return rcpp_result_gen;
END_RCPP
}
// timeseries1D
NumericVector timeseries1D(const unsigned int& N, const double& startpoint, const double& d13, const double& d12, const double& d11, const double& d10, const double& d22, const double& d21, const double& d20, const double& sf, double dt);
RcppExport SEXP _Langevin_timeseries1D(SEXP NSEXP, SEXP startpointSEXP, SEXP d13SEXP, SEXP d12SEXP, SEXP d11SEXP, SEXP d10SEXP, SEXP d22SEXP, SEXP d21SEXP, SEXP d20SEXP, SEXP sfSEXP, SEXP dtSEXP) {
......@@ -84,6 +103,7 @@ END_RCPP
static const R_CallMethodDef CallEntries[] = {
{"_Langevin_Langevin1D", (DL_FUNC) &_Langevin_Langevin1D, 6},
{"_Langevin_Langevin2D", (DL_FUNC) &_Langevin_Langevin2D, 6},
{"_Langevin_kernel1D", (DL_FUNC) &_Langevin_kernel1D, 4},
{"_Langevin_timeseries1D", (DL_FUNC) &_Langevin_timeseries1D, 11},
{"_Langevin_timeseries2D", (DL_FUNC) &_Langevin_timeseries2D, 11},
{NULL, NULL, 0}
......
// Copyright (C) 2021 Philip Rinn
//
// 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>.
// [[Rcpp::depends(RcppArmadillo)]]
#define ARMA_NO_DEBUG
#include <RcppArmadillo.h>
using namespace Rcpp;
// Epanechnikov kernel
arma::vec kappa(const arma::vec& data) {
arma::vec res = arma::zeros<arma::vec>(data.n_elem);
arma::uvec hit = arma::find(data%data <= 5);
res(hit) = 0.75 * (1 - data(hit)%data(hit)/5) / std::sqrt(5);
return res;
}
// Nadaraya-Watson estimator
// [[Rcpp::export(".kernel1D")]]
List kernel1D(const arma::vec& data, const int& bins, const double& sf, const double& h) {
arma::vec ests = arma::linspace<arma::vec>(data.min(), data.max(), bins);
arma::vec ka_sum(bins);
arma::vec wa_sum(bins);
arma::vec wa2_sum(bins);
arma::vec wa4_sum(bins);
arma::vec phi(bins);
ka_sum.fill(NA_REAL);
wa_sum.fill(NA_REAL);
wa2_sum.fill(NA_REAL);
wa4_sum.fill(NA_REAL);
phi.fill(NA_REAL);
arma::vec data_min = data.head(data.n_elem - 1);
arma::vec data_inc = data.tail(data.n_elem - 1) - data_min;
for(int i = 0; i < bins; i++) {
arma::vec data_kappa = kappa((ests(i) - data_min)/h)/h;
ka_sum(i) = arma::sum(data_kappa);
wa_sum(i) = arma::sum(data_inc%data_kappa);
wa2_sum(i) = arma::sum(data_inc%data_inc%data_kappa);
wa4_sum(i) = arma::sum(data_inc%data_inc%data_inc%data_inc%data_kappa);
}
arma::vec M1 = wa_sum / ka_sum;
arma::vec M2 = wa2_sum / ka_sum;
arma::vec M4 = wa4_sum / ka_sum;
arma::vec eD1 = sf * arma::sqrt((M2 - arma::square(M1)) / ka_sum);
arma::vec eD2 = (sf/2) * arma::sqrt((M4 - arma::square(M2)) / ka_sum);
arma::vec D1 = sf * M1;
arma::vec D2 = (sf/2) * M2;
arma::vec D4 = (sf/24) * M4;
List ret;
ret["D1"] = D1;
ret["eD1"] = eD1;
ret["D2"] = D2;
ret["eD2"] = eD2;
ret["D4"] = D4;
ret["mean_bin"] = ests;
ret.attr("class") = CharacterVector::create("Langevin");
return ret;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment