# R/Kernel_Silverman.R In distr6: The Complete R6 Probability Distributions Interface

#' @title Silverman Kernel
#'
#' @description Mathematical and statistical functions for the Silverman kernel defined by the pdf,
#' \deqn{f(x) = exp(-|x|/\sqrt{2})/2 * sin(|x|/\sqrt{2} + \pi/4)}
#' over the support \eqn{x \in R}{x \epsilon R}.
#'
#' @template param_decorators
#' @template class_kernel
#' @template class_distribution
#' @template method_pdfsquared2Norm
#'
#' @details The cdf and quantile functions are omitted as no closed form analytic expressions could
#' be found, decorate with FunctionImputation for numeric results.
#'
#' @export
Silverman <- R6Class("Silverman",
inherit = Kernel,
lock_objects = FALSE,
public = list(
name = "Silverman",
short_name = "Silv",
description = "Silverman Kernel",

#' @description
#' Creates a new instance of this [R6][R6::R6Class] class.
initialize = function(decorators = NULL) {
super$initialize( decorators = decorators, support = Reals$new()
)
},

#' @description
#' The squared 2-norm of the pdf is defined by
#' \deqn{\int_a^b (f_X(u))^2 du}
#' where X is the Distribution, \eqn{f_X} is its pdf and \eqn{a, b}
#' are the distribution support limits.
pdfSquared2Norm = function(x = 0, upper = Inf) {

xl <- length(x)
ul <- length(upper)
len <- max(xl, ul)

ret <- numeric(len)
for (i in seq(len)) {

xi <- x[ifelse(i %% xl == 0, xl, i %% xl)]
ui <- upper[ifelse(i %% ul == 0, ul, i %% ul)]

if (xi >= 0) {
if (ui == Inf) {
ret[i] <- (exp(-xi / sqrt(2)) * (3 * sin(xi / sqrt(2)) + 3 *
cos(xi / sqrt(2)))) / 2^ (7 / 2) +
(xi * exp(-xi / sqrt(2)) * sin(xi / sqrt(2))) / 8
} else if (ui <= 0) {
ret[i] <- (exp((2 * ui - xi) / sqrt(2)) *
(2 * cos(xi / sqrt(2)) - sin((2 * ui - xi) / sqrt(2)) +
cos((2 * ui - xi) / sqrt(2)))) / 2^ (9 / 2)
} else if (ui >= 0 & ui <= xi) {
ret[i] <- (exp(-xi / sqrt(2)) * ((2 * ui + sqrt(2)) * sin(xi / sqrt(2)) +
sqrt(2) * sin((2 * ui - xi) / sqrt(2)))) / 16 +
(exp(- xi / sqrt(2)) * (sin(xi / sqrt(2)) + 3 * cos(xi / sqrt(2)))) / 2^ (9 / 2)
} else if (ui >= xi) {
ret[i] <- (exp(- xi  / sqrt(2)) / (8 * sqrt(2))) *
((3 + sqrt(2) * xi) * sin(xi / sqrt(2)) + 3 * cos(xi / sqrt(2))) +
(exp((xi - 2 * ui) / sqrt(2)) / (16 * sqrt(2))) *
(sin((xi - 2 * ui) / sqrt(2)) -
cos((xi - 2 * ui) / sqrt(2)) - 2 * cos(xi / sqrt(2)))
}
} else if (xi <= 0) {
if (ui == Inf) {
ret[i] <- (exp(xi / sqrt(2)) * (-3 * sin(xi / sqrt(2)) + 3 *
cos(xi / sqrt(2)))) / 2^ (7 / 2) +
(xi * exp(xi / sqrt(2)) * sin(xi / sqrt(2))) / 8
} else if (ui <= xi) {
ret[i] <- (exp((2 * ui - xi) / sqrt(2)) / (16 * sqrt(2))) *
(sin((xi - 2 * ui) / sqrt(2)) + cos((xi - 2 * ui) / sqrt(2)) +
2 * cos(xi / sqrt(2)))
} else if (ui >= xi & ui <= 0) {
ret[i] <- (- (exp(xi / sqrt(2)) * (2 * sin((xi - 2 * ui) / sqrt(2)) -
(2^ (3 / 2) * (xi - ui) - 3) * sin(xi / sqrt(2)) -
3 * cos(xi / sqrt(2))))) / 2^ (9 / 2)
} else if (ui >= 0) {
ret[i] <- (exp((xi - 2 * ui) / sqrt(2)) *
(sin((xi - 2 * ui) / sqrt(2)) -
cos((xi - 2 * ui) / sqrt(2)) +
exp(sqrt(2) * ui) *
(2^ (3 / 2) * xi - 6) * sin(xi / sqrt(2)) +
2 * (3 * exp(sqrt(2) * ui) - 1) * cos(xi / sqrt(2)))) / 2^ (9 / 2)
}
}
}
return(ret)
},

#' @description
#' The squared 2-norm of the cdf is defined by
#' \deqn{\int_a^b (F_X(u))^2 du}
#' where X is the Distribution, \eqn{F_X} is its pdf and \eqn{a, b}
#' are the distribution support limits.
cdfSquared2Norm = function(x = 0, upper = 0) {

xl <- length(x)
ul <- length(upper)
len <- max(xl, ul)

ret <- numeric(len)
for (i in seq(len)) {

xi <- x[ifelse(i %% xl == 0, xl, i %% xl)]
ui <- upper[ifelse(i %% ul == 0, ul, i %% ul)]

if (xi >= 0) {
if (ui >= -Inf & ui <= 0) {
ret[i] <- (exp((2 * ui - xi) / sqrt(2)) *
(2 * cos(xi / sqrt(2)) + sin((2 * ui - xi) / sqrt(2)) +
cos((2 * ui - xi) / sqrt(2)))) / (2^ (9 / 2))
} else if (ui >= 0 & ui <= xi) {
ret[i] <- (exp(- xi / sqrt(2)) *
(sqrt(2) * (sin((xi - 2 * ui) / sqrt(2)) +
3 * sin(xi / sqrt(2))) -
2 * (ui + 2^ (3 / 2)) * cos(xi / sqrt(2)) +
2^ (5 / 2) * exp(ui / sqrt(2)) *
(sin((ui - xi) / sqrt(2)) +
cos((ui - xi) / sqrt(2))))) / 16 +
(exp(- xi / sqrt(2)) * (3 * cos(xi / sqrt(2)) -
sin(xi / sqrt(2)))) / (2^ (9 / 2))
} else if (ui >= xi) {
ret[i] <- (8 * exp((xi - ui) / (sqrt(2))) *
(sin((xi - ui) / (sqrt(2))) +
cos((xi - ui) / (sqrt(2)))) -
exp((xi - 2 * ui) / (sqrt(2))) *
(sin((xi - 2 * ui) / (sqrt(2))) +
cos((xi - 2 * ui) / (sqrt(2))) +
2 * cos((xi) / (sqrt(2))))) / (16 * sqrt(2)) -
(8 * exp((-ui) / (sqrt(2))) * (sin((ui) / (sqrt(2))) -
cos((ui) / (sqrt(2))))) / (16 * sqrt(2)) +
(exp((- xi) / (sqrt(2))) * (10 * sin((xi) / (sqrt(2))) -
(10 + 2 * sqrt(2) * xi) * cos((xi) / (sqrt(2)))) -
16 * sqrt(2) * (xi - ui)) / (16 * sqrt(2))
}
} else if (xi <= 0) {
if (ui <= xi) {
ret[i] <-  (exp((2 * ui - xi) / sqrt(2)) * (2 * cos(xi / sqrt(2)) +
sin((2 * ui - xi) / sqrt(2)) +
cos((2 * ui - xi) / sqrt(2)))) / (2^ (9 / 2))
} else if (ui >= xi  & ui <= 0) {
ret[i] <- (exp(xi / sqrt(2)) *
(sqrt(2) * (sin((xi - 2 * ui) / sqrt(2)) -
3 * sin(xi / sqrt(2))) -
2 * (- xi + ui + 2^ (3 / 2)) *
cos(xi / sqrt(2))) +
2^ (5 / 2) * exp(ui / sqrt(2)) * (sin(ui / sqrt(2)) +
cos(ui / sqrt(2)))) / 16 +
(exp(xi / sqrt(2)) * (sin(xi / sqrt(2)) +
3 * cos(xi / sqrt(2)))) / (2^ (9 / 2))
} else if (ui >= 0) {
ret[i] <- (exp(- ui / sqrt(2)) *
(exp(ui / sqrt(2)) *
(256 * exp((xi - ui) / sqrt(2)) *
(sin((xi - ui) / sqrt(2)) +
cos((xi - ui) / sqrt(2))) -
32 * exp((xi - 2 * ui) / sqrt(2)) *
sin((xi - 2 * ui) / sqrt(2)) -
32 * exp((xi - 2 * ui) / sqrt(2)) *
cos((xi - 2 * ui) / sqrt(2)) -
320 * exp(xi / sqrt(2)) * sin(xi / sqrt(2)) -
64 * exp((xi - 2 * ui) / sqrt(2)) * cos(xi / sqrt(2)) +
2^ (13 / 2) * xi * exp(xi /  sqrt(2)) * cos(xi / sqrt(2)) -
320 * exp(xi / sqrt(2)) * cos(xi / sqrt(2)) +
2^ (19 / 2) * ui) - 256 * sin(ui / sqrt(2)) +
256 * cos(ui / sqrt(2)))) / 2 ^ (19 / 2)
}
}
}
return(ret)
},

#' @description
#' The variance of a distribution is defined by the formula
#' \deqn{var_X = E[X^2] - E[X]^2}
#' where \eqn{E_X} is the expectation of distribution X. If the distribution is multivariate the
#' covariance matrix is returned.
#' @param ... Unused.
variance = function(...) {
return(0)
}
),

private = list(
.isQuantile = 0L,
.pdf = function(x, log = FALSE) {
C_SilvermanKernelPdf(x, log)
},
.cdf = function(x, lower.tail = TRUE, log.p = FALSE) {
C_SilvermanKernelCdf(x, lower.tail, log.p)
}
)
)
.distr6$kernels <- rbind( .distr6$kernels,
data.table::data.table(
ShortName = "Silv", ClassName = "Silverman",
Support = "\u211D", Packages = "-"
)
)


## Try the distr6 package in your browser

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

distr6 documentation built on March 28, 2022, 1:05 a.m.