24.10.22.1

make sure `na_action` is also effective for `method = "MVKE"`

DESCRIPTION

@@ -50,5 +50,5 @@ Suggests:
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.0
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr

NEWS.md

@@ -4,6 +4,7 @@
 - Debug: 
 	- Changed the function form of `find_eqs()` according to the new setting of the `MVKE()` function; added `linear_interp` to `sim_vf_options`.
 	- Fixed a typo in the `fit_2d_vf()` function.
+	- The parameter `na_action` in `fit_2d_vf()` was not effective for `method = "MVKE"` in the previous version. Now it is fixed. The `vector_position` parameter is now also effective for `method = "MVKE"`.
 
 # fitlandr 0.1.0
 

R/2d_landscape.R

@@ -9,6 +9,7 @@
 #' @param method The method used to estimate the gradient. Currently only "MVKE" is supported.
 #' @param ... Additional arguments passed to [MVKE()].
 #' @inheritParams stats::integrate
+#' @inheritParams fit_2d_vf
 #' @return A `2d_MVKE_landscape` object, which contains the following components:
 #' \itemize{
 #'   \item `dist`: A data frame containing the estimated potential landscape. The data frame has two columns: `x` and `U`, where `x` is the position and `U` is the potential.
@@ -25,10 +26,58 @@
 #' summary(l)
 #' plot(l)
 #'
+#' # different behaviors for different `na_action` choices
 #'
-fit_2d_ld <- function(data, x, lims, n = 200L, method = c("MVKE"), subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL, ...) {
+#' l1 <- fit_2d_ld(data.frame(x = c(1,2,1,2,NA,NA,NA,10,11,10,11)), "x")
+#' plot(l1)
+#'
+#' l2 <- fit_2d_ld(data.frame(x = c(1,2,1,2,NA,NA,NA,10,11,10,11)), "x", na_action = "omit_vectors")
+#' plot(l2)
+#'
+#'
+fit_2d_ld <- function(data, x, lims, n = 200L, vector_position = "start", na_action = "omit_data_points",
+											method = c("MVKE"), subdivisions = 100L, rel.tol = .Machine$double.eps^0.25, abs.tol = rel.tol, stop.on.error = TRUE, keep.xy = FALSE, aux = NULL, ...) {
+	d <- data
+	# extract useful data for construction
+	if (is.data.frame(d)) {
+		d_raw <- d[, c(x), drop = FALSE] %>% as.matrix()
+	} else if (is.matrix(d)) {
+		d_raw <- d[, c(x), drop = FALSE]
+	} else {
+		rlang::abort("`d` must be a data frame or a matrix.")
+	}
+
+	if (na_action != "omit_data_points" & na_action != "omit_vectors") {
+		rlang::abort('`na_action` must be either "omit_data_points" or "omit_vectors".')
+	}
+
+	if (na_action == "omit_data_points" & any(is.na(d_raw))) {
+		d_raw <- stats::na.omit(d_raw)
+		rlang::inform("NA(s) found in the data. Those data points were omitted.")
+	}
+
+	if (vector_position == "start") {
+		x_mat <- d_raw[1:(nrow(d_raw) - 1), , drop = FALSE]
+	} else if (vector_position == "middle") {
+		x_mat <- d_raw[1:(nrow(d_raw) - 1), , drop = FALSE] + 0.5 * v_mat
+	} else if (vector_position == "end") {
+		x_mat <- d_raw[2:nrow(d_raw), , drop = FALSE]
+	} else {
+		rlang::abort('`vector_position` must be one of "start", "middle", or "end".')
+	}
+
+	v_mat <- diff(d_raw)
+
+	data_vectors <- cbind(x_mat, v_mat) %>%
+		`colnames<-`(c("x", "vx"))
+
+	if (any(is.na(data_vectors))) {{ if (na_action == "omit_vectors") {
+		data_vectors <- stats::na.omit(data_vectors)
+		rlang::inform("NA(s) found in the data. Those vectors were omitted.")
+	} }}
+
 	lims <- determine_lims(data, x, lims)
-  MVKEresult <- MVKE(data[, x, drop = FALSE], ...)
+  MVKEresult <- MVKE(data_vectors[,1, drop = FALSE], data_vectors[,2, drop = FALSE], ...)
 
   xseq <- seq(lims[1], lims[2], length.out = n)
   Useq <- vector("numeric", length = n)

R/fit_vectorfield.R

@@ -7,7 +7,7 @@
 #' @param x,y Characters to indicate the name of the two variables.
 #' @param lims The limits of the range for the vector field estimation as `c(<xl>, <xu>, <yl>, <yu>)`. If missing, the range of the data extended by 10% for both sides will be used.
 #' @param n The number of equally spaced points in each axis, at which the vectors are to be estimated.
-#' @param vector_position Only useful if `method == "VFC"`. One of "start", "middle", or "end", representing the position of the vectors. If "start", for example, the starting point of a vector is regarded as the position of the vector.
+#' @param vector_position One of "start", "middle", or "end", representing the position of the vectors. If "start", for example, the starting point of a vector is regarded as the position of the vector.
 #' @param na_action One of "omit_data_points" or "omit_vectors". If using "omit_data_points", then only the `NA` points are omitted, and the points before and after an `NA` will form a vector. If using "omit_vectors", then the vectors will be omitted if either of its points is `NA`.
 #' @param method One of "MVKE" or "VFC".
 #' @param ... Other parameters to be passed to [MVKE()] or [SparseVFC::SparseVFC()].
@@ -80,7 +80,7 @@ fit_2d_vf <- function(data, x, y,
   if (method == "VFC") {
     VFCresult <- SparseVFC::SparseVFC(original_vectors_normalized[, 1:2], original_vectors_normalized[, 3:4], ...)
   } else if (method == "MVKE") {
-    MVKEresult <- MVKE(original_vectors_normalized[, 1:2], ...)
+    MVKEresult <- MVKE(original_vectors_normalized[, 1:2], original_vectors_normalized[, 3:4], ...)
   }
 
   lims <- determine_lims(d_raw, c(x, y), lims)

R/kernal_estimator.R

@@ -3,21 +3,33 @@
 #' See references for details.
 #'
 #' @param d The dataset. Should be a matrix or a data frame, with each row representing a random vector.
+#' @param v The vectors corresponding to the dataset. Should be a matrix or a data frame with the same shape as `d`. If missing, then the vectors will be calculated from the dataset.
 #' @param h The bandwidth for the kernel estimator.
 #' @param kernel The type of kernel estimator used. "exp" by default ([exp()]), and if "Gaussian" then [stats::dnorm()] will be used.
 #'
 #' @return A function(x), which then returns the \eqn{\mu} and \eqn{a} estimators at the position \eqn{x}.
 #' @references Bandi, F. M., & Moloche, G. (2018). On the functional estimation of multivariate diffusion processes. Econometric Theory, 34(4), 896-946. https://doi.org/10.1017/S0266466617000305
 #' @export
-MVKE <- function(d, h = 0.2, kernel = c("exp", "Gaussian")) {
+MVKE <- function(d, v, h = 0.2, kernel = c("exp", "Gaussian")) {
   if (is.data.frame(d)) d <- as.matrix(d)
   if (!is.matrix(d)) stop("`d` should be a data.frame or a matrix.")
+  if (any(is.na(d))) stop("There are missing values in `d`.")
+  if (missing(v)) {
+    v <- diff(d)
+    d <- d[1:(nrow(d) - 1), , drop = FALSE]
+  } else {
+    if (is.data.frame(v)) v <- as.matrix(v)
+    if (!is.matrix(v)) stop("`v` should be a data.frame or a matrix.")
+    if (any(is.na(v))) stop("There are missing values in `v`.")
+    if (!all(dim(v) == dim(d))) stop("`v` should have the same shape as `d`.")
+  }
+
 
-  d <- stats::na.omit(d)
+  # d <- stats::na.omit(d)
   dim <- ncol(d)
 
-  temp_d <- d[1:(nrow(d) - 1), , drop = FALSE]
-  temp_diff <- diff(d)
+  temp_d <- d
+  temp_diff <- v
   temp_norm <- apply(temp_diff, MARGIN = 1, FUN = function(x) norm(x, "2"))
   temp_diff_tcrossprod <- apply(temp_diff,
     MARGIN = 1,