R/DNAmf.R

In DNAmf: Diffusion Non-Additive Model with Tunable Precision

#' @title Fitting a Diffusion Non-Additive model for multi-fidelity computer experiments with tuning parameters
#'
#' @description The function fits DNA models for multi-fidelity computer experiments with tuning parameters.
#' Available kernel choices include nonseparable squared exponential kernel, and nonseparable Matern kernel with smoothness parameter 1.5 and 2.5.
#' The function returns a fitted model object of class \code{DNAmf}, produced by \code{DNAmf_internal}.
#'
#' @seealso \code{\link{predict.DNAmf}} for prediction.
#'
#' @details The \code{DNAmf} function internally calls \code{DNAmf_internal} to fit the DNA model with nonseparable kernel.
#'
#' The model structure is:
#' \eqn{\begin{cases}
#' & f_1(\bm{x}) = W_1(\bm{x}),\\
#' & f_l(\bm{x}) = W(t_l, \bm{x}, f_{l-1}(\bm{x})),
#' \end{cases}}
#' where \eqn{W(t, \bm{x}, y) \sim GP(\alpha, \tau^2 K((t, \bm{x}, y), (t', \bm{x}', y')))} is a GP model.
#' Hyperparameters \eqn{(\alpha, \tau^2, \bm{\theta})} are estimated by
#' maximizing the log-likelihood via an optimization algorithm "L-BFGS-B".
#' For \code{constant=FALSE}, \eqn{\alpha=0}.
#'
#' The nonseparable covariance kernel is defined as:
#' \deqn{K((t, \bm{x}, y), (t', \bm{x}', y'))=
#' \left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ - \left(\frac{\beta(d+1)}{2}+\delta \right) }
#' \prod^d_{j=1}\phi(x_j,x'_j;\theta_{j})\phi(y,y';\theta_{y}),}
#' where \eqn{\phi(\cdot, \cdot)} depens on the chosen kernel:
#'
#' \itemize{
#' \item For nonseparable squared exponential kernel(\code{kernel = "sqex"}):
#' \deqn{\phi(x, x';\theta) = \exp \left( -\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{-\beta}
#' \frac{ (x-x')^2 }{\theta} \right)}
#'
#' \item For nonseparable Matern kernel with smoothness parameter of \eqn{\nu=1.5} (\code{kernel = "matern1.5"}):
#' \deqn{\phi(x,x';\theta) = \left( 1+\frac{1}{\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ \frac{\beta}{2} }}\frac{\sqrt{3}|x- x'|}{\theta} \right)
#' \exp \left( -\frac{1}{\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ \frac{\beta}{2} }}\frac{\sqrt{3}|x- x'|}{\theta} \right)}
#'
#' \item For nonseparable Matern kernel with smoothness parameter of \eqn{\nu=2.5} (\code{kernel = "matern2.5"}):
#' \deqn{\phi(x, x';\theta) = \left( 1+\frac{1}{\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ \frac{\beta}{2} }}\frac{\sqrt{5}|x- x'|}{\theta}+
#' \frac{1}{3}\left(\frac{1}{\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ \frac{\beta}{2} }}\frac{\sqrt{5}|x- x'|}{\theta} \right)^2 \right) }
#' \deqn{\times \exp \left( -\frac{1}{\left(\frac{(t-t')^2}{\theta_t} + 1\right)^{ \frac{\beta}{2} }}\frac{\sqrt{5}|x- x'|}{\theta} \right)}
#' }
#'
#' When the input locations are not nested, the internal \code{makenested} function constructs nested designs as
#' \eqn{\mathcal{X}^*_L = \mathcal{X}_L} and
#' \eqn{\mathcal{X}^*_l = \mathcal{X}_l \cup \mathcal{X}^*_{l+1}} for \eqn{l = 1, \dots, L-1}.
#' The function \code{\link{imputer}} then imputes pseudo outputs \eqn{\widetilde{\mathbf{y}}_l := f_l(\widetilde{\mathcal{X}}_l)}
#' at pseudo inputs \eqn{\widetilde{\mathcal{X}}_l := \mathcal{X}^*_l \setminus \mathcal{X}_l},
#' using a stochastic EM algorithm.
#'
#' For further details, see Heo, Boutelet, and Sung (2025+, <arXiv:2506.08328>).
#'
#' @param X A list of input locations for all fidelity levels \eqn{1,\ldots,L} combined.
#' @param y A list of response values for all fidelity levels \eqn{1,\ldots,L} combined.
#' @param kernel A character specifying the kernel type to be used. Choices are \code{"sqex"}(squared exponential), \code{"matern1.5"}, or \code{"matern2.5"}. Default is \code{"sqex"}.
#' @param t A vector of tuning parameters for each fidelity level.
#' @param constant A logical indicating for constant mean of GP (\code{constant=TRUE}) or zero mean (\code{constant=FALSE}). Default is \code{TRUE}.
#' @param init Optional vector of initial parameter values for optimization. Default is \code{NULL}.
#' @param n.iter Number of iterations for the stochastic EM algorithm for non-nested designs. Default is 50.
#' @param multi.start Number of random starting points for optimization. Default is 10.
#' @param g Nugget term for numerical stability. Default is \code{sqrt(.Machine$double.eps)}.
#' @param burn.ratio Fraction of iterations to discard as burn-in. Default is 0.75.
#' @param ... Additional arguments for compatibility with \code{DNAmf_internal}.
#'
#' @return A fitted model object of class \code{DNAmf}.
#'
#' @usage DNAmf(X, y, kernel = "sqex", t, constant = TRUE, init=NULL,
#' n.iter=50, multi.start=10, g = sqrt(.Machine$double.eps), burn.ratio = 0.75, ...)
#' @export
#' @examples
#' ### Non-Additive example ###
#' library(RNAmf)
#'
#' ### Non-Additive Function ###
#' fl <- function(x, t){
#'   term1 <- sin(10 * pi * x / (5+t))
#'   term2 <- 0.2 * sin(8 * pi * x)
#'   term1 + term2
#' }
#'
#' ### training data ###
#' n1 <- 13; n2 <- 10; n3 <- 7; n4 <- 4; n5 <- 1;
#' m1 <- 2.5; m2 <- 2.0; m3 <- 1.5; m4 <- 1.0; m5 <- 0.5;
#' d <- 1
#'
#' ### fix seed to reproduce the result ###
#' set.seed(1)
#'
#' ### generate initial nested design ###
#' NestDesign <- NestedX(c(n1,n2,n3,n4,n5),d)
#'
#' X1 <- NestDesign[[1]]
#' X2 <- NestDesign[[2]]
#' X3 <- NestDesign[[3]]
#' X4 <- NestDesign[[4]]
#' X5 <- NestDesign[[5]]
#'
#' y1 <- fl(X1, t=m1)
#' y2 <- fl(X2, t=m2)
#' y3 <- fl(X3, t=m3)
#' y4 <- fl(X4, t=m4)
#' y5 <- fl(X5, t=m5)
#'
#' ### fit a DNAmf ###
#' fit.DNAmf <- DNAmf(X=list(X1, X2, X3, X4, X5), y=list(y1, y2, y3, y4, y5), kernel="sqex",
#'                    t=c(m1,m2,m3,m4,m5), multi.start=10, constant=TRUE)
#'

DNAmf <- function(X, y, kernel="sqex", t, constant = TRUE, init=NULL, n.iter=50, multi.start=10,
                  g = sqrt(.Machine$double.eps), burn.ratio = 0.75, ...) {

  g <- sqrt(.Machine$double.eps)
  X1 <- X[[1]]; y1 <- y[[1]]
  nn <- unlist(lapply(X, nrow))
  X <- do.call(rbind, X[-1])
  y <- do.call(rbind, y[-1])
  lvl <- rep(seq_along(nn[-1]), times = nn[-1])
  idxs  <- split(seq_len(nrow(X)), lvl)
  X_list <- lapply(idxs, function(i) X[i, , drop = FALSE])
  y_list <- lapply(idxs, function(i) y[i, , drop = FALSE])

  # check whether the design is nested #
  nested <- all(unlist(lapply(X_list, checknested, XX1=X1)))

  # make designs nested
  if(!nested) XX <- makenested(c(list(X1), X_list)) else XX <- X
  model <- DNAmf_internal(X1, y1, XX, c(list(y1), y_list), kernel=kernel, t=t, nn=nn, nested = nested, constant = constant,
                          init=init, n.iter = n.iter, multi.start = multi.start, trace = TRUE, g = g, burn.ratio = 0.75, ...)
  return(model)
}