R/RcppExports.R

In ldsr: Linear Dynamical System Reconstruction

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Implement Kalman smoothing
#'
#' Estimate the hidden state and expected log-likelihood given the observations, exogeneous input and system parameters. This is an internal function and should not be called directly.
#'
#' @param y Observation matrix (may need to be normalized and centered before hand) (q rows, T columns)
#' @param u Input matrix for the state equation (m_u rows, T columns)
#' @param v Input matrix for the output equation (m_v rows, T columns)
#' @param theta A list of system parameters (A, B, C, D, Q, R)'
#' @param stdlik Boolean, whether the likelihood is divided by the number of observations. Standardizing the likelihood this way may speed up convergence in the case of long time series.
#' @return A list of fitted elements (X, Y, V, J, and lik)
#' @section Note: This code only works on one dimensional state and output at the moment. Therefore, transposing is skipped, and matrix inversion is treated as /, and log(det(Sigma)) is treated as log(Sigma).
Kalman_smoother <- function(y, u, v, theta, stdlik = TRUE) {
    .Call(`_ldsr_Kalman_smoother`, y, u, v, theta, stdlik)
}

#' Maximizing expected likelihood using analytical solution
#'
#' @inheritParams Kalman_smoother
#' @param fit result of [Kalman_smoother]
#' @return An object of class `theta`: a list of
Mstep <- function(y, u, v, fit) {
    .Call(`_ldsr_Mstep`, y, u, v, fit)
}

#' Learn LDS model
#'
#' Estimate the hidden state and model parameters given observations and exogenous inputs using the EM algorithm. This is the key backend routine of this package.
#'
#' @inheritParams Kalman_smoother
#' @param theta0 A vector of initial values for the parameters
#' @param niter Maximum number of iterations, default 1000
#' @param tol Tolerance for likelihood convergence, default 1e-5. Note that the log-likelihood is normalized
#' @return A list of model results
#' * theta: model parameters (A, B, C, D, Q, R, mu1, V1) resulted from Mstep
#' * fit: results of Estep
#' * liks : vector of loglikelihood over the iteration steps
#' @section Note: This code only works on one dimensional state and output at the moment. Therefore, transposing is skipped, and matrix inversion is treated as /, and log(det(Sigma)) is treated as log(Sigma).
LDS_EM <- function(y, u, v, theta0, niter = 1000L, tol = 1e-5) {
    .Call(`_ldsr_LDS_EM`, y, u, v, theta0, niter, tol)
}

#' State propagation
#'
#' This function propagates the state trajectory based on the exogenous inputs only
#' (without measurement update), and calculates the corresponding log-likelihood
#'
#' @param theta A list of system parameters (A, B, C, D, Q, R)'
#' @param u Input matrix for the state equation (m_u rows, T columns)
#' @param v Input matrix for the output equation (m_v rows, T columns)
#' @param y Observations
#' @param stdlik Boolean, whether the likelihood is divided by the number of observations. Standardizing the likelihood this way may speed up convergence in the case of long time series.
#' @section Note: This code only works on one dimensional state and output at the moment. Therefore, transposing is skipped, and matrix inversion is treated as /, and log(det(Sigma)) is treated as log(Sigma).
#' @return A list of predictions and log-likelihood (X, Y, V, lik)
propagate <- function(theta, u, v, y, stdlik = TRUE) {
    .Call(`_ldsr_propagate`, theta, u, v, y, stdlik)
}

#' Nash-Sutcliffe Efficiency
#'
#' @param yhat Model outputs
#' @param y Observations
#' @return NSE
#' @examples
#' NSE(rnorm(100), rnorm(100))
#' @export
NSE <- function(yhat, y) {
    .Call(`_ldsr_NSE`, yhat, y)
}

#' Normalized root-mean-square error
#'
#' RMSE is normalized by the normalization constant
#' @param yhat Model outputs
#' @param y Observations
#' @param normConst The normalization constant
#' @return normalized RMSE
#' @examples
#' x <- rnorm(100)
#' y <- rnorm(100)
#' nRMSE(x, y, sd(y))
#' @export
nRMSE <- function(yhat, y, normConst) {
    .Call(`_ldsr_nRMSE`, yhat, y, normConst)
}

#' Pearson's correlation
#'
#' Calculate the Pearson's correlation using the numerically stable formulation (see References). Internal function.
#' @param x First variable
#' @param y Second variable
#' @return Pearson's correlation
#' @section Reference John D. Cook's article at https://www.johndcook.com/blog/2008/11/05/how-to-calculate-pearson-correlation-accurately/
corr <- function(x, y) {
    .Call(`_ldsr_corr`, x, y)
}

#' Kling-Gupta Efficiency
#'
#' @param yhat Model outputs
#' @param y Observations
#' @return KGE
#' @examples
#' KGE(rnorm(100), rnorm(100))
#' @export
KGE <- function(yhat, y) {
    .Call(`_ldsr_KGE`, yhat, y)
}

#' Reduction of Error
#'
#' @param yhat Model outputs in the validation set
#' @param y Observations in the validation set
#' @param yc_bar Mean observations in the calibration set
#' @return RE
#' @examples
#' x <- rnorm(100)
#' y <- rnorm(100)
#' yc_bar <- mean(x[1:50])
#' RE(x[51:100], y[51:100], yc_bar)
#' @export
RE <- function(yhat, y, yc_bar) {
    .Call(`_ldsr_RE`, yhat, y, yc_bar)
}