R/bernstein_polynomials.R In beyondWhittle: Bayesian Spectral Inference for Time Series

Documented in betaBasis_kcoarsened_bernsteincoarsened_bernstein_idbList

```#' Construct coarsened Bernstein polynomial basis of degree l on omega
#' @param omega numeric vector in [0,1] of evaluation points
#' @param l positive integer for the degree
#' @details See Appendix E.1 in Ghosal/Van der Vaart, Fundamentals, (2017)
#' @keywords internal
coarsened_bernstein <- function(omega, l) {
res <- matrix(NA, nrow=l, ncol=length(omega))
for (i in 1:l) {
res[i,] <- coarsened_bernstein_i(omega, l, i)
}
res
}

#' Helping function for \code{coarsened_bernstein}
#' @keywords internal
coarsened_bernstein_i <- function(omega, l, i) {
k <- l^2
b_tmp <- 0 * omega
for (j in ((i-1)*l+1):(i*l)) {
b_tmp <- b_tmp + dbeta(omega, j, k+1-j)
}
b_tmp <- b_tmp / l
b_tmp
}

#' Construct Bernstein polynomial basis of degree k on omega
#' @param omega numeric vector in [0,1] of evaluation points
#' @param k positive integer for the degree
#' @param coarsened bool flag indicating whether coarsened or standard Bernstein polynomials are used
#' @keywords internal
betaBasis_k <- function(omega, k, coarsened) {
if (coarsened) {
basis <- coarsened_bernstein(omega, k)
} else {
N <- length(omega)
basis <- matrix(dbeta(omega,
rep(1:k, each = N),
rep(k:1, each = N)),
ncol = N,
byrow = TRUE)
}
basis
}

#' Construct Bernstein polynomial basises of degree up to kmax on omega
#' @param n positive integer determining the number of the (equidistant) evaluation points in [0,1]
#' @param kmax positive integer for the largest degree
#' @param bernstein_l,bernstein_r left and right truncation
#' @param coarsened bool flag indicating whether coarsened or standard Bernstein polynomials are used
#' @param verbose debugging parameter
#' @return A list of length kmax, where the k-th list element is a matrix containing the polynomial basis of degree k
#' @keywords internal
dbList <- function(n, kmax, bernstein_l=0, bernstein_r=1, coarsened=F, verbose=F) {
db.list <- vector("list", kmax)
omega <- omegaFreq(n); NN <- length(omega)
omega_for_dblist <- seq(bernstein_l, bernstein_r, length.out=NN)
if (verbose) {
if (coarsened) {
cat("Using coarsened Bernstein polynomials on (", bernstein_l , ",", bernstein_r, ")\n")
} else {
cat("Using standard Bernstein polynomials on (", bernstein_l , ",", bernstein_r, ")\n")
}
}
if (verbose) cat("Precomputing polynomial basis functions ")
for (kk in 1:kmax) {
db.list[[kk]] <- betaBasis_k(omega_for_dblist, kk, coarsened)
if (verbose && !(kk%%10)) cat(".")
}
if (verbose) cat(" done!\n")
return(db.list)
}
```

Try the beyondWhittle package in your browser

Any scripts or data that you put into this service are public.

beyondWhittle documentation built on May 31, 2023, 6:51 p.m.