# R/utils.R In hermiter: Efficient Sequential and Batch Estimation of Univariate and Bivariate Probability Density Functions and Cumulative Distribution Functions along with Quantiles (Univariate) and Spearman's Correlation (Bivariate)

#' Convenience function to output Hermite normalization factors
#'
#' The method returns numeric normalization factors that, when multiplied by
#' the physicist Hermite polynomials \eqn{H_k(x)}, yield orthonormal
#' Hermite functions \eqn{h_k(x)} for \eqn{k=0,\dots,N}.
#'
#' @author Michael Stephanou <michael.stephanou@gmail.com>
#'
#' @param N An integer number.
#' @return A numeric vector of length N+1
#' @export
hermite_normalization_N <- function(N){
if (N < 0) {
stop("N must be >= 0.")
}
return(hermite_normalization(N))
}

#' Convenience function to output physicist Hermite polynomials
#'
#'
#' The method calculates the physicist version of Hermite polynomials,
#' \eqn{H_k(x)} from \eqn{k=0,\dots,N} for the vector of values, x.
#'
#' @param N An integer number.
#' @param x A numeric vector.
#' @return A numeric matrix with N+1 rows and length(x) columns.
#' @export
hermite_polynomial_N <- function(N,x){
if (N < 0) {
stop("N must be >= 0.")
}
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (length(x)<1) {
stop("x must contain at least one value.")
}
return(hermite_polynomial(N,x))
}

#' Convenience function to output orthonormal Hermite functions
#'
#'
#' The method calculates the orthonormal Hermite functions, \eqn{h_k(x)}
#' from \eqn{k=0,\dots,N} for the vector of values, x.
#'
#' @param N An integer number.
#' @param x A numeric vector.
#' @param normalization_hermite A numeric vector. A vector of normalization
#' values necessary in the calculation of the Hermite functions.
#' @return A numeric matrix with N+1 rows and length(x) columns.
#' @export
hermite_function_N <- function(N,x, normalization_hermite=NULL){
if (N < 0) {
stop("N must be >= 0.")
}
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (length(x)<1) {
stop("x must contain at least one value.")
}
if (is.null(normalization_hermite)){
normalization_hermite <- hermite_normalization(N)
} else {
if (length(normalization_hermite) != (N+1)){
stop("Normalization vector must be of length N+1.")
}
}
return(hermite_function(N,x,normalization_hermite))
}

#' Convenience function to output the sum of orthonormal Hermite functions
#'
#'
#' The method calculates the sum of orthonormal Hermite functions,
#' \eqn{\sum_{i} h_k(x_{i})} from \eqn{k=0,\dots,N} for the vector of values,
#' x.
#'
#'
#' @param N An integer number.
#' @param x A numeric vector.
#' @param normalization_hermite A numeric vector of normalization factors
#' generated by the hermite_normalization function.
#' @return A numeric vector of length N+1.
#' @export
hermite_function_sum_N <- function(N,x, normalization_hermite=NULL){
if (N < 0) {
stop("N must be >= 0.")
}
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (length(x)<1) {
stop("x must contain at least one value.")
}
if (is.null(normalization_hermite)){
normalization_hermite <- hermite_normalization(N)
} else {
if (length(normalization_hermite) != (N+1)){
stop("Normalization vector must be of length N+1.")
}
}
return(hermite_function_sum(N,x,normalization_hermite))
}

#' Convenience function to output a definite integral of the orthonormal
#' Hermite functions
#'
#' The method calculates \eqn{\int_{-\infty}^{x} h_k(t) dt}
#' for \eqn{k=0,\dots,N} and the vector of values x.
#'
#' @param N An integer number.
#' @param x A numeric vector.
#' @param hermite_function_matrix A numeric matrix. A matrix of Hermite
#' function values.
#' @param normalization_hermite A numeric vector. A vector of normalization
#' values necessary in the calculation of the Hermite functions.
#' @return A numeric matrix with N+1 rows and length(x) columns.
#' @export
hermite_int_lower <- function(N,x,hermite_function_matrix=NULL,
normalization_hermite=NULL){
if (N < 0) {
stop("N must be >= 0.")
}
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (length(x)<1) {
stop("x must contain at least one value.")
}
if (is.null(hermite_function_matrix)){
if (is.null(normalization_hermite)){
normalization_hermite <- hermite_normalization(N)
} else {
if (length(normalization_hermite) != (N+1)){
stop("Normalization vector must be of length N+1.")
}
}
hermite_function_matrix <- hermite_function_N(N,x,normalization_hermite)
} else {
if (nrow(hermite_function_matrix) != (N+1)){
stop("Hermite function matrix must have N+1 rows.")
}
}
return(hermite_integral_val(N,x,hermite_function_matrix))
}

#' Convenience function to output a definite integral of the orthonormal
#' Hermite functions
#'
#' The method calculates \eqn{\int_{x}^{\infty} h_k(t) dt}
#' for \eqn{k=0,\dots,N} and the vector of values x.
#'
#' @param N An integer number.
#' @param x A numeric vector.
#' @param hermite_function_matrix A numeric matrix. A matrix of Hermite
#' function values.
#' @param normalization_hermite A numeric vector. A vector of normalization
#' values necessary in the calculation of the Hermite functions.
#' @return A numeric matrix with N+1 rows and length(x) columns.
#' @export
hermite_int_upper <- function(N,x, hermite_function_matrix=NULL,
normalization_hermite=NULL){
if (N < 0) {
stop("N must be >= 0.")
}
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (length(x)<1) {
stop("x must contain at least one value.")
}
if (is.null(hermite_function_matrix)){
if (is.null(normalization_hermite)){
normalization_hermite <- hermite_normalization(N)
} else {
if (length(normalization_hermite) != (N+1)){
stop("Normalization vector must be of length N+1.")
}
}
hermite_function_matrix <- hermite_function_N(N,x,normalization_hermite)
} else {
if (nrow(hermite_function_matrix) != (N+1)){
stop("Hermite function matrix must have N+1 rows.")
}
}
return(hermite_integral_val_upper(N,x,hermite_function_matrix))
}

#' Convenience function to output the integral of the orthonormal Hermite
#' functions on the full domain
#'
#' The method calculates \eqn{\int_{-\infty}^{\infty} h_k(t) dt}
#' for \eqn{k=0,\dots,N}.
#'
#' @param N An integer number.
#' @return A numeric matrix with N+1 rows and 1 columns.
#' @export
hermite_int_full <- function(N){
if (N < 0) {
stop("N must be >= 0.")
}
return(hermite_int_full_domain(N))
}

#' Calculates \eqn{\int_{-\infty}^{\infty} f(x) e^{-x^2} dx} using
#' Gauss-Hermite quadrature with 100 terms.
#'
#' @param f A function.
#' @return A numeric value.
#' @export
result <- sum(f(root_x_serialized)*weight_w_serialized)
return(result)
}

# Helper method for merging univariate hermite_estimator objects.
integrand_coeff_univar <- function(t,hermite_est_current,
hermite_estimator_merged, current_k,
dimension = NA){
normalization_hermite <- hermite_est_current$normalization_hermite_vec t <- sqrt(2) * t herm_mod <- hermite_polynomial_N(hermite_est_current$N_param, t) *
normalization_hermite
if (is.na(dimension)){
original_sd <- sqrt(hermite_est_current$running_variance / (hermite_est_current$num_obs-1))
original_mean <- hermite_est_current$running_mean new_sd <- sqrt(hermite_estimator_merged$running_variance /
(hermite_estimator_merged$num_obs-1)) new_mean <- hermite_estimator_merged$running_mean
original_coeff_vec <- hermite_est_current$coeff_vec } else { if (dimension==1){ original_sd <- sqrt(hermite_est_current$running_variance_x /
(hermite_est_current$num_obs-1)) original_mean <- hermite_est_current$running_mean_x
new_sd <- sqrt(hermite_estimator_merged$running_variance_x / (hermite_estimator_merged$num_obs-1))
new_mean <- hermite_estimator_merged$running_mean_x original_coeff_vec <- hermite_est_current$coeff_vec_x
} else if (dimension==2) {
original_sd <- sqrt(hermite_est_current$running_variance_y / (hermite_est_current$num_obs-1))
original_mean <- hermite_est_current$running_mean_y new_sd <- sqrt(hermite_estimator_merged$running_variance_y /
(hermite_estimator_merged$num_obs-1)) new_mean <- hermite_estimator_merged$running_mean_y
original_coeff_vec <- hermite_est_current\$coeff_vec_y
}
}
return(sqrt(2) * hermite_function_N(current_k-1,((t*original_sd +
original_mean) -  new_mean)/new_sd,
normalization_hermite[1:current_k])[current_k, ] *
as.vector(crossprod(herm_mod, original_coeff_vec)))
}

# Helper method for estimating Hermite series sums, with or without series
# acceleration. These series acceleration techniques are drawn from the below
# reference:
#
# Boyd, John P., and Dennis W. Moore. "Summability methods for
# Hermite functions." Dynamics of atmospheres and oceans 10.1 (1986): 51-62.
#
# Note that h_input is a numeric matrix of N+1 rows and length(x) columns,
# coeffs is a numeric vector of length N+1
#
# This helper method is intended for internal use by the
# hermite_estimator_univar class.
series_calculate <- function(h_input, coeffs, accelerate_series = TRUE){
if (!is.numeric(coeffs)) {
stop("coeffs must be numeric.")
}
if (!is.numeric(h_input)) {
stop("h_input must be numeric.")
}
if (is.null(nrow(h_input))) {
stop("h_input must be numeric matrix.")
}
if (nrow(h_input) != length(coeffs)) {
stop("Number of rows of h_input and length of coeffs must match.")
}
sz <- nrow(h_input) - 1
if (length(coeffs) < 3 | accelerate_series == FALSE){
return(as.numeric(crossprod(h_input,coeffs)))
}
if (length(coeffs) >=3 & length(coeffs) < 6){
result <- h_input[1,] * coeffs[1] + 1/2 * h_input[2,] * coeffs[2]
return(as.numeric(result))
}
if (length(coeffs) >=6 & length(coeffs) < 12){
result <- crossprod(h_input[1:(sz-3),,drop=FALSE], coeffs[1:(sz-3)])
result <- result + 15/16*h_input[sz-2,]*coeffs[sz-2]+
11/16*h_input[sz-1,]*coeffs[sz-1]+
5/16*h_input[sz,]*coeffs[sz]+
1/16*h_input[sz+1,]*coeffs[sz+1]
return(as.numeric(result))
}
if (length(coeffs) >= 12){
result <- crossprod(h_input[1:(sz-8),,drop=FALSE], coeffs[1:(sz-8)])
result <- result + 255/256*h_input[sz-7,]*coeffs[sz-7]+
247/256*h_input[sz-6,]*coeffs[sz-6]+
219/256*h_input[sz-5,]*coeffs[sz-5]+
163/256*h_input[sz-4,]*coeffs[sz-4]+
93/256*h_input[sz-3,]*coeffs[sz-3]+
37/256*h_input[sz-2,]*coeffs[sz-2]+
9/256*h_input[sz-1,]*coeffs[sz-1]+
1/256*h_input[sz+1,]*coeffs[sz+1]
return(as.numeric(result))
}
}


## Try the hermiter package in your browser

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

hermiter documentation built on Nov. 17, 2021, 1:07 a.m.