Nothing
R/kedd.R
In shazam: Immunoglobulin Somatic Hypermutation Analysis
Defines functions
kernel_fun_conv
kernel_fun_der
A2_kM
A3_kMr
h.ucv
# This file contains functions adapted from the CRAN package kedd,
# scheduled for archival on 2022-05-25
# Original code: https://github.com/cran/kedd
# Modification: Simplify to use only kernel="gaussian"
#
# kedd: Kernel Estimator and Bandwidth Selection for Density and Its Derivatives
# Smoothing techniques and computing bandwidth selectors of the nth derivative of a probability density for one-dimensional data.
#
# Version: 1.0.3
# Depends: R (≥ 2.15.0)
# Published: 2015-10-31
# Author: Arsalane Chouaib Guidoum
# Maintainer: Arsalane Chouaib Guidoum <acguidoum at usthb.dz>
# License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
# NeedsCompilation: no
# Classification/MSC: 62G05, 62G07, 65D10, 68N15
h.ucv <- function(x,deriv.order=0,lower=0.1*hos,upper=2*hos,tol=0.1 * lower,
kernel="gaussian",...)
{
if (!is.numeric(x) || length(dim(x)) >=1 || length(x) < 3L)
stop("argument 'x' must be numeric and need at least 3 data points")
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
if (missing(kernel)) kernel <- "gaussian"
if (kernel != "gaussian") {stop("Expecting gaussian kernel.")}
name <- deparse(substitute(x))
x <- x[!is.na(x)]
x <- sort(x)
n <- length(x)
hos <- ((243 *(2*r+1)*A3_kMr(kernel,r))/(35* A2_kM(kernel)^2))^(1/(2*r+5)) * sd(x,na.rm = TRUE) * n^(-1/(2*r+5))
if (!is.numeric(upper)){
stop("argument 'upper' must be numeric. Default 2*hos (Oversmoothing) boundary was used")
upper= 2*hos
}
if (!is.numeric(lower)){
stop("argument 'lower' must be numeric. Default 0.1*hos boundary was used")
lower=0.1*hos
}
if (lower < 0 | lower >= upper){
stop("the boundaries must be positive and 'lower' must be smaller than 'upper'. Default boundaries were used")
upper=2*hos
lower=0.1*hos
}
R_Kr1 <- A3_kMr(kernel,r)
fucv <- function(h)
{
D <- kernel_fun_der(kernel, outer(x,x,"-")/h,deriv.order=2*r)
diag(D) <- 0
D <- ((-1)^r / ((n-1)*h^(2*r+1)))* colSums(D)
D1 <- mean(D)
D2 <- kernel_fun_conv(kernel,outer(x,x,"-")/h,deriv.order=r)
diag(D2) <- 0
D3 <- ((-1)^r / ((n-1)*h^(2*r+1)))* colSums(D2)
D4 <- mean(D3)
(1/(n*h^(2*r+1)))* R_Kr1 + D4 - 2*D1
}
obj <- optimize(fucv , c(lower, upper),tol=tol)
structure(list(x=x, data.name=name,n=n, kernel=kernel, deriv.order=r, h = obj$minimum ,
min.ucv=obj$objective),class="h.ucv")
}
A3_kMr <- function(kernel,r) {
xKr <- integrate(function(x) kernel_fun_der(kernel,x,deriv.order=r)^2, -Inf, Inf)$value
return(xKr)
}
A2_kM <-function(kernel) {
xKr <- 1
return(xKr)
}
####
#### r(th) derivative of Kernel functions K^r(x)
kernel_fun_der <- function(kernel,u,deriv.order=0)
{
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
Kr <- expression( dnorm(X) )
if (r == 0) {
DKr <- Kr
K <- function(X) eval(DKr);fx <- K(u)
} else {
if (r == 1){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)/sqrt(pi)}
else if (r == 2){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^2-1)/sqrt(pi) }
else if (r == 3){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^2-3)/sqrt(pi)}
else if (r == 4){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^4-6*u^2+3)/sqrt(pi)}
else if (r == 5){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^4-10*u^2+15)/sqrt(pi)}
else if (r == 6){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^6-15*u^4+45*u^2-15)/sqrt(pi)}
else if (r == 7){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^6-21*u^4+105*u^2-105)/sqrt(pi)}
else if (r == 8){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^8-28*u^6+210*u^4-420*u^2+105)/sqrt(pi)}
else if (r == 9){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^8-36*u^6+378*u^4-1260*u^2+945)/sqrt(pi)}
else if (r == 10){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^10-45*u^8+630*u^6-3150*u^4+4725*u^2-945)/sqrt(pi)}
else if (r == 11){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^10-55*u^8+990*u^6-6930*u^4+17325*u^2-10395)/sqrt(pi)}
else if (r == 12){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^12-66*u^10+1485*u^8-13860*u^6+51975*u^4-62370*u^2+10395)/sqrt(pi)}
else if (r == 13){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^12-78*u^10+2145*u^8-25740*u^6+135135*u^4-270270*u^2+135135)/sqrt(pi)}
else if (r == 14){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^14-91*u^12+3003*u^10-45045*u^8+315315*u^6-945945*u^4+945945*u^2-135135)/sqrt(pi)}
else if (r == 15){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^14-105*u^12+4095*u^10-75075*u^8+675675*u^6-2837835*u^4+4729725*u^2-2027025)/sqrt(pi)}
else if (r == 16){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^16-120*u^14+5460*u^12-120120*u^10+1351350*u^8-7567560*u^6+18918900*u^4-16216200*u^2+2027025)/sqrt(pi)}
else if (r == 17){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^16-136*u^14+7140*u^12-185640*u^10+2552550*u^8-18378360*u^6+64324260*u^4-91891800*u^2+34459425)/sqrt(pi)}
else if (r == 18){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^18-153*u^16+9180*u^14-278460*u^12+4594590*u^10-41351310*u^8+192972780*u^6-413513100*u^4+310134825*u^2-34459425)/sqrt(pi)}
else if (r == 19){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^18-171*u^16+11628*u^14-406980*u^12+7936110*u^10-87297210*u^8+523783260*u^6-1571349780*u^4+1964187225*u^2-654729075)/sqrt(pi)}
else if (r == 20){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^20-190*u^18+14535*u^16-581400*u^14+13226850*u^12-174594420*u^10+1309458150*u^8-5237832600*u^6+9820936125*u^4-6547290750*u^2+654729075)/sqrt(pi)}
else if (r == 21){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^20-210*u^18+17955*u^16-813960*u^14+21366450*u^12-333316620*u^10+3055402350*u^8-15713497800*u^6+41247931725*u^4-45831035250*u^2+13749310575)/sqrt(pi)}
else if (r == 22){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^22-231*u^20+21945*u^18-1119195*u^16+33575850*u^14-611080470*u^12+6721885170*u^10-43212118950*u^8+151242416325*u^6-252070693875*u^4+151242416325*u^2-13749310575)/sqrt(pi)}
else if (r == 23){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^22-253*u^20+26565*u^18-1514205*u^16+51482970*u^14-1081142370*u^12+14054850810*u^10-110430970650*u^8+496939367925*u^6-1159525191825*u^4+1159525191825*u^2-316234143225)/sqrt(pi)}
else if (r == 24){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^24-276*u^22+31878*u^20-2018940*u^18+77224455*u^16-1853386920*u^14+28109701620*u^12-265034329560*u^10+1490818103775*u^8-4638100767300*u^6+6957151150950*u^4-3794809718700*u^2+316234143225)/sqrt(pi)}
else if (r == 25){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^24-300*u^22+37950*u^20-2656500*u^18+113565375*u^16-3088978200*u^14+54057118500*u^12-602350749000*u^10+4141161399375*u^8-16564645597500*u^6+34785755754750*u^4-31623414322500*u^2+7905853580625)/sqrt(pi)}
else if (r == 26){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^26-325*u^24+44850*u^22-3453450*u^20+164038875*u^18-5019589575*u^16+100391791500*u^14-1305093289500*u^12+10767019638375*u^10-53835098191875*u^8+150738274937250*u^6-205552193096250*u^4+102776096548125*u^2-7905853580625)/sqrt(pi)}
else if (r == 27){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^26-351*u^24+52650*u^22-4440150*u^20+233107875*u^18-7972289325*u^16+180705224700*u^14-2710578370500*u^12+26428139112375*u^10-161505294575625*u^8+581419060472250*u^6-1109981842719750*u^4+924984868933125*u^2-213458046676875)/sqrt(pi)}
else if (r == 28){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^28-378*u^26+61425*u^24-5651100*u^22+326351025*u^20-12401338950*u^18+316234143225*u^16-5421156741000*u^14+61665657928875*u^12-452214824811750*u^10+2034966711652875*u^8-5179915266025500*u^6+6474894082531875*u^4-2988412653476250*u^2+213458046676875)/sqrt(pi)}
else if (r == 29){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^28-406*u^26+71253*u^24-7125300*u^22+450675225*u^20-18928359450*u^18+539458244325*u^16-10480903032600*u^14+137561852302875*u^12-1192202719958250*u^10+6557114959770375*u^8-21459648959248500*u^6+37554385678684875*u^4-28887988983603750*u^2+6190283353629375)/sqrt(pi)}
else if (r == 30){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^30-435*u^28+82215*u^26-8906625*u^24+614557125*u^22-28392539175*u^20+899097073875*u^18-19651693186125*u^16+294775397791875*u^14-2980506799895625*u^12+19671344879311125*u^10-80473683597181875*u^8+187771928393424375*u^6-216659917377028125*u^4+92854250304440625*u^2-6190283353629375)/sqrt(pi)}
else {K <- function(X) (-1)^r * Hermite(X,r) * eval(Kr);fx <- K(u)}
}
return(fx)
}
#### Hermite Polynomial
Hermite <-function (x, n, prob = TRUE) {
if (any(n < 0 || n != round(n))) stop("Argument 'n' must be a vector of non-negative integers")
if ((length(n) != 1) && (length(x) != length(n)) && (length(x) != 1))
stop(paste("Argument 'n' must be either a vector of same length",
"as argument 'x',\n a single integer or 'x' must be a ",
"single value!", sep = ""))
H <- function(x, n) {
if (n <= 1) {
return(switch(n + 1, 1, x))
}
else {
return(x * Recall(x, n - 1) - (n - 1) * Recall(x, n - 2))
}
}
scale <- 1
if (!prob) {
x <- sqrt(2) * x
scale <- 2^(n/2)
}
scale * mapply(H, x, n)
}
#### Kernels density convolution
kernel_fun_conv <- function(kernel,u,deriv.order=0)
{
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
if (r==0) {fx <- dnorm(u,mean=0,sd=sqrt(2))}
else if (r==1) {fx <- (1/8)*exp(-(1/4)*u^2)*(u^2-2)/sqrt(pi)}
else if (r==2) {fx <- (1/32)*exp(-(1/4)*u^2)*(12-12*u^2+u^4)/sqrt(pi)}
else if (r==3) {fx <- (1/128)*exp(-(1/4)*u^2)*(u^6-30*u^4+180*u^2-120)/sqrt(pi)}
else if (r==4) {fx <- (1/512)*exp(-(1/4)*u^2)*(u^8-56*u^6+840*u^4-3360*u^2+1680)/sqrt(pi)}
else if (r==5) {fx <- (1/2048)*exp(-(1/4)*u^2)*(u^10-90*u^8+2520*u^6-25200*u^4+75600*u^2-30240)/sqrt(pi)}
else if (r==6) {fx <- (1/8192)*exp(-(1/4)*u^2)*(u^12-132*u^10+5940*u^8-110880*u^6+831600*u^4-1995840*u^2+665280)/sqrt(pi)}
else if (r==7) {fx <- (1/32768)*exp(-(1/4)*u^2)*(u^14-182*u^12+12012*u^10-360360*u^8+5045040*u^6-30270240*u^4+60540480*u^2-17297280)/sqrt(pi)}
else if (r==8) {fx <- (1/131072)*exp(-(1/4)*u^2)*(u^16-240*u^14+21840*u^12-960960*u^10+21621600*u^8-242161920*u^6+1210809600*u^4-2075673600*u^2+518918400)/sqrt(pi)}
else if (r==9) {fx <- (1/524288)*exp(-(1/4)*u^2)*(u^18-306*u^16+36720*u^14-2227680*u^12+73513440*u^10-1323241920*u^8+12350257920*u^6-52929676800*u^4+79394515200*u^2-17643225600)/sqrt(pi)}
else if (r==10) {fx <- (1/2097152)*exp(-(1/4)*u^2)*(u^20-380*u^18+58140*u^16-4651200*u^14+211629600*u^12-5587021440*u^10+83805321600*u^8-670442572800*u^6+2514159648000*u^4-3352212864000*u^2+670442572800)/sqrt(pi)}
else if (r==11) {fx <- (1/8388608)*exp(-(1/4)*u^2)*(u^22-462*u^20+87780*u^18-8953560*u^16+537213600*u^14-19554575040*u^12+430200650880*u^10-5531151225600*u^8+38718058579200*u^6-129060195264000*u^4+154872234316800*u^2-28158588057600)/sqrt(pi)}
else if (r==12) {fx <- (1/33554432)*exp(-(1/4)*u^2)*(u^24-552*u^22+127512*u^20-16151520*u^18+1235591280*u^16-59308381440*u^14+1799020903680*u^12-33924394183680*u^10+381649434566400*u^8-2374707592857600*u^6+7124122778572800*u^4-7771770303897600*u^2+1295295050649600)/sqrt(pi)}
else if (r==13) {fx <- (1/134217728)*exp(-(1/4)*u^2)*(u^26-650*u^24+179400*u^22-27627600*u^20+2624622000*u^18-160626866400*u^16+6425074656000*u^14-167051941056000*u^12+2756357027424000*u^10-27563570274240000*u^8+154355993535744000*u^6-420970891461120000*u^4+420970891461120000*u^2-64764752532480000)/sqrt(pi)}
else if (r==14) {fx <- (1/536870912)*exp(-(1/4)*u^2)*(u^28-756*u^26+245700*u^24-45208800*u^22+5221616400*u^20-396842846400*u^18+20238985166400*u^16-693908062848000*u^14+15786408429792000*u^12-231533990303616000*u^10+2083805912732544000*u^8-10608466464820224000*u^6+26521166162050560000*u^4-24481076457277440000*u^2+3497296636753920000)/sqrt(pi)}
else if (r==15) {fx <- (1/2147483648)*exp(-(1/4)*u^2)*(u^30-870*u^28+328860*u^26-71253000*u^24+9832914000*u^22-908561253600*u^20+57542212728000*u^18-2515416727824000*u^16+75462501834720000*u^14-1526019481546560000*u^12+20143457156414592000*u^10-164810104007028480000*u^8+769113818699466240000*u^6-1774878043152614400000*u^4+1521324036987955200000*u^2-202843204931727360000)/sqrt(pi)}
else if (r==16) {fx <- (1/8589934592)*exp(-(1/4)*u^2)*(u^32-992*u^30+431520*u^28-108743040*u^26+17670744000*u^24-1950850137600*u^22+150215460595200*u^20-8154553575168000*u^18+311911674250176000*u^16-8317644646671360000*u^14+151381132569418752000*u^12-1816573590833025024000*u^10+13624301931247687680000*u^8-58689300626913116160000*u^6+125762787057670963200000*u^4-100610229646136770560000*u^2+12576278705767096320000)/sqrt(pi)}
else if (r==17) {fx <- (1/34359738368)*exp(-(1/4)*u^2)*(u^34-1122*u^32+556512*u^30-161388480*u^28+30502422720*u^26-3965314953600*u^24+364808975731200*u^22-24077392398259200*u^20+1143676138917312000*u^18-38884988723188608000*u^16+933239729356526592000*u^14-15440875522080712704000*u^12+169849630742887839744000*u^10-1175882058989223505920000*u^8+4703528235956894023680000*u^6-9407056471913788047360000*u^4+7055292353935341035520000*u^2-830034394580628357120000)/sqrt(pi)}
else if (r==18) {fx <- (1/137438953472)*exp(-(1/4)*u^2)*(u^36-1260*u^34+706860*u^32-233735040*u^30+50837371200*u^28-7686610525440*u^26+832716140256000*u^24-65665615631616000*u^22+3792189302725824000*u^20-160114659448423680000*u^18+4899508579121764608000*u^16-106898368999020318720000*u^14+1621291929818474833920000*u^12-16462348825849129082880000*u^10+105829385309030115532800000*u^8-395096371820379097989120000*u^6+740805697163210808729600000*u^4-522921668585795864985600000*u^2+58102407620643984998400000)/sqrt(pi)}
else if (r==19) {fx <- (1/549755813888)*exp(-(1/4)*u^2)*(u^38-1406*u^36+885780*u^34-331281720*u^32+82157866560*u^30-14295468781440*u^28+1801229066461440*u^26-167256984742848000*u^24+11540731947256512000*u^22-592424239959167616000*u^20+22512121118448369408000*u^18-626246278385927367168000*u^16+12524925567718547343360000*u^14-175348957948059662807040000*u^12+1653290174938848249323520000*u^10-9919741049633089495941120000*u^8+34719093673715813235793920000*u^6-61268988835969082180812800000*u^4+40845992557312721453875200000*u^2-4299578163927654889881600000)/sqrt(pi)}
else if (r==20) {fx <- (1/2199023255552)*exp(-(1/4)*u^2)*(u^40-1560*u^38+1096680*u^36-460605600*u^34+129199870800*u^32-25633254366720*u^30+3716821883174400*u^28-401416763382835200*u^26+32615112024855360000*u^24-2000393537524462080000*u^22+92418181433630148096000*u^20-3192628085889041479680000*u^18+81412016190170557731840000*u^16-1502991068126225681203200000*u^14+19538883885640933855641600000*u^12-171942178193640217929646080000*u^10+967174752339226225854259200000*u^8-3185987419470392273402265600000*u^6+5309979032450653789003776000000*u^4-3353670967863570814107648000000*u^2+335367096786357081410764800000)/sqrt(pi)}
else if (r==21) {fx <- (1/8796093022208)*exp(-(1/4)*u^2)*(u^42-1722*u^40+1343160*u^38-629494320*u^36+198290710800*u^34-44496435503520*u^32+7356744003248640*u^30-914338183260902400*u^28+86404958318155276800*u^26-6240358100755658880000*u^24+344467767161712370176000*u^22-14467646220791919547392000*u^20+458142130325077452334080000*u^18-10783960913805669262632960000*u^16+184867901379525758787993600000*u^14-2243063870071579206627655680000*u^12+18505276928090528454678159360000*u^10-97969113148714562407119667200000*u^8+304792796462667527488816742400000*u^6-481251783888422411824447488000000*u^4+288751070333053447094668492800000*u^2-27500101936481280675682713600000)/sqrt(pi)}
else if (r==22) {fx <- (1/35184372088832)*exp(-(1/4)*u^2)*(u^44-1892*u^42+1629012*u^40-847086240*u^38+297750813360*u^36-75033204966720*u^34+14031209328776640*u^32-1988422807735203840*u^30+216240980341203417600*u^28-18164242348661087078400*u^26+1180675752662970660096000*u^24-59248455951814527670272000*u^22+2281065554144859315305472000*u^20-66677300813465118447390720000*u^18+1457375289208594731778682880000*u^16-23318004627337515708458926080000*u^14+265242302635964241183720284160000*u^12-2059528467526310578603004559360000*u^10+10297642337631552893015022796800000*u^8-30350945837229840105728488243200000*u^6+45526418755844760158592732364800000*u^4-26015096431911291519195847065600000*u^2+2365008766537390138108713369600000)/sqrt(pi)}
else if (r==23) {fx <- (1/140737488355328)*exp(-(1/4)*u^2)*(u^46-2070*u^44+1958220*u^42-1124018280*u^40+438367129200*u^38-123268836731040*u^36+25886455713518400*u^34-4149229044366806400*u^32+514504401501483993600*u^30-49735425478476786048000*u^28+3759998166172845025228800*u^26-222181709819304478763520000*u^24+10220358651688006023121920000*u^22-363215822852296829437102080000*u^20+9858715191705199656149913600000*u^18-201117789910786072985458237440000*u^16+3016766848661791094781873561600000*u^14-32297150968026234073547116953600000*u^12+236845773765525716539345524326400000*u^10-1121901033626174446765320904704000000*u^8+3141322894153288450942898533171200000*u^6-4487604134504697787061283618816000000*u^4+2447784073366198792942518337536000000*u^2-212850788988365112429784203264000000)/sqrt(pi)}
else if (r==24) {fx <- (1/562949953421312)*exp(-(1/4)*u^2)*(u^48-2256*u^46+2334960*u^44-1472581440*u^42+633946309920*u^40-197791248695040*u^38+46349082610871040*u^36-8342834869956787200*u^34+1170082590511439404800*u^32-128969103309705321062400*u^30+11220311987944362932428800*u^28-771141442080539852446924800*u^26+41770161446029242007541760000*u^24-1773625316785241660627927040000*u^22+58529635453912974800721592320000*u^20-1482750764832462028284947005440000*u^18+28357608377420836290949611479040000*u^16-400342706504764747636935691468800000*u^14+4047909587992621337217905324851200000*u^12-28122319242896106132250710677913600000*u^10+126550436593032477595128198050611200000*u^8-337467830914753273587008528134963200000*u^6+460183405792845373073193447456768000000*u^4-240095689978875846820796581281792000000*u^2+20007974164906320568399715106816000000)/sqrt(pi)}
else if (r==25) {fx <- (1/2251799813685248)*exp(-(1/4)*u^2)*(u^50-2450*u^48+2763600*u^46-1906884000*u^44+901956132000*u^42-310633691860800*u^40+80764759883808000*u^38-16222178913804864000*u^36+2554993178924266080000*u^34-318522482972558504640000*u^32+31597430310877803660288000*u^30-2499069488223971744040960000*u^28+157441377758110219874580480000*u^26-7872068887905510993729024000000*u^24+310384430437417290609887232000000*u^22-9559840457472452550784526745600000*u^20+227046210864970748081132510208000000*u^18-4086831795569473465460385183744000000*u^16+54491090607592979539471802449920000000*u^14-521967288977995909272835160309760000000*u^12+3444984107254773001200712058044416000000*u^10-14764217602520455719431623105904640000000*u^8+37581644806415705467644131542302720000000*u^6-49019536704020485392579302011699200000000*u^4+24509768352010242696289651005849600000000*u^2-1960781468160819415703172080467968000000)/sqrt(pi)}
else if (r==26) {fx <- (1/9007199254740992)*exp(-(1/4)*u^2)*(u^52-2652*u^50+3248700*u^48-2443022400*u^46+1264264092000*u^44-478397532412800*u^42+137300091802473600*u^40-30598306173122688000*u^38+5377652309926312416000*u^36-752871323389683738240000*u^34+84472162484322515430528000*u^32-7617853198586175937007616000*u^30+552294356897497755433052160000*u^28-32118041062654484854414417920000*u^26+1491194763623243939669240832000000*u^24-54875967301335376979828062617600000*u^22+1584543555826059010292535308083200000*u^20-35419208894935436700656671592448000000*u^18+602126551213902423911163417071616000000*u^16-7605809067965083249404169478799360000000*u^14+69212862518482257569577942257074176000000*u^12-435052278687602761865918494187323392000000*u^10+1779759321903829480360575658039050240000000*u^8-4333327044635410908704010297834209280000000*u^6+5416658805794263635880012872292761600000000*u^4-2599996226781246545222406178700525568000000*u^2+199999709752403580401723552207732736000000)/sqrt(pi)}
else if (r==27) {fx <- (1/36028797018963968)*exp(-(1/4)*u^2)*(u^54-2862*u^52+3795012*u^50-3099259800*u^48+1747982527200*u^46-723664766260800*u^44+228195622960905600*u^42-56136123248382777600*u^40+10946544033434641632000*u^38-1710093434556567348288000*u^36+215471772754127485884288000*u^34-21978120820921003560197376000*u^32+1816857987862802960976316416000*u^30-121589726880049121234568867840000*u^28+6565845251522652546666718863360000*u^26-284519960899314943688891150745600000*u^24+9815938651026365557266744700723200000*u^22-266762568045540052203366826572595200000*u^20+5631654214294734435404410783199232000000*u^18-90699273135483617749144721034682368000000*u^16+1088391277625803412989736652416188416000000*u^14-9432724406090296245911050987606966272000000*u^12+56596346436541777475466305925641797632000000*u^10-221463964316902607512694240578598338560000000*u^8+516749250072772750862953228016729456640000000*u^6-620099100087327301035543873620075347968000000*u^4+286199584655689523554866403209265545216000000*u^2-21199969233754779522582696534019670016000000)/sqrt(pi)}
else if (r==28) {fx <- (1/144115188075855872)*exp(-(1/4)*u^2)*(u^56-3080*u^54+4407480*u^52-3896212320*u^50+2386430046000*u^48-1076757236755200*u^46+371481246680544000*u^44-100406074102798464000*u^42+21612407450627369376000*u^40-3746150624775410691840000*u^38+526708777843422743272704000*u^36-60332096371155696047600640000*u^34+5641051010703057580450659840000*u^32-430455584816725624600542658560000*u^30+26749739913610806671605150924800000*u^28-1348186891645984656248899606609920000*u^26+54770092473118126660111546518528000000*u^24-1778417120303600348022445510483968000000*u^22+45646039421125742265909434769088512000000*u^20-912920788422514845318188695381770240000000*u^18+13967688062864477133368287039341084672000000*u^16-159630720718451167238494709021040967680000000*u^14+1320581416852641474427547138264975278080000000*u^12-7578989001067333679323314010912032030720000000*u^10+28421208754002501297462427540920120115200000000*u^8-63663507608965602906315837691661069058048000000*u^6+73457893394960311045749043490378156605440000000*u^4-32647952619982360464777352662390291824640000000*u^2+2331996615713025747484096618742163701760000000)/sqrt(pi)}
else if (r==29) {fx <- (1/576460752303423488)*exp(-(1/4)*u^2)*(u^58-3306*u^56+5091240*u^54-4857042960*u^52+3220219482480*u^50-1577907546415200*u^48+593293237452115200*u^46-175445285932268352000*u^44+41492810122981465248000*u^42-7938957670197120350784000*u^40+1238477396550750774722304000*u^38-158299929050032326296323584000*u^36+16621492550253394261113976320000*u^34-1434562664721869873920760110080000*u^32+101649011671721065352099573514240000*u^30-5895642676959821790421775263825920000*u^28+278569116486351579597428881215774720000*u^26-10651172100948736866960516046485504000000*u^24+326635944429094597253455825425555456000000*u^22-7942410859275879154268241649821401088000000*u^20+150905806326241703931096591346606620672000000*u^18-2198913177896664828710264616764839329792000000*u^16+23988143758872707222293795819252792688640000000*u^14-189819224526731857150324819091478620405760000000*u^12+1044005734897025214326786505003132412231680000000*u^10-3758420645629290771576431418011276684034048000000*u^8+8095059852124626277241544592639672857919488000000*u^6-8994510946805140308046160658488525397688320000000*u^4+3854790405773631560591211710780796599009280000000*u^2-265847614191284935213187014536606662000640000000)/sqrt(pi)}
else if (r==30) {fx <- (1/2305843009213693952)*exp(-(1/4)*u^2)*(u^60-3540*u^58+5851620*u^56-6007663200*u^54+4298483019600*u^52-2279915393595840*u^50+930965452384968000*u^48-300036865797212544000*u^46+77634539025028745760000*u^44-16320505315039376330880000*u^42+2810391015249780604177536000*u^40-398564543980877976592450560000*u^38+46698479069759536257415457280000*u^36-4526160279069001206487959705600000*u^34+362739416651101382405677913548800000*u^32-23989166754526171423095499349360640000*u^30+1304410942277360571130817777121484800000*u^28-58007921903628505398523425853167206400000*u^26+2094730513186584917168901489142149120000000*u^24-60857433856789203909328085368761384960000000*u^22+1405806722091830610305478772018387992576000000*u^20-25438407352137887234099139684142258913280000000*u^18+353825120443372431528833488333978692157440000000*u^16-3692088213322147111605219008702386352947200000000*u^14+27998335617692948929672910815993096509849600000000*u^12-147831212061418770348672969108443549572005888000000*u^10+511723426366449589668483354606150748518481920000000*u^8-1061352291723006556349446957701645996927221760000000*u^6+1137163169703221310374407454680334996707737600000000*u^4-470550277118574335327341015729793791741132800000000*u^2+31370018474571622355156067715319586116075520000000)/sqrt(pi)}
else if (r>=31) {fx <- NA}
return(fx)
}
Try the shazam package in your browser
Any scripts or data that you put into this service are public.
shazam documentation built on Oct. 3, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions kernel_fun_conv kernel_fun_der A2_kM A3_kMr h.ucv
# This file contains functions adapted from the CRAN package kedd,
# scheduled for archival on 2022-05-25
# Original code: https://github.com/cran/kedd
# Modification: Simplify to use only kernel="gaussian"
#
# kedd: Kernel Estimator and Bandwidth Selection for Density and Its Derivatives
# Smoothing techniques and computing bandwidth selectors of the nth derivative of a probability density for one-dimensional data.
#
# Version: 1.0.3
# Depends: R (≥ 2.15.0)
# Published: 2015-10-31
# Author: Arsalane Chouaib Guidoum
# Maintainer: Arsalane Chouaib Guidoum <acguidoum at usthb.dz>
# License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
# NeedsCompilation: no
# Classification/MSC: 62G05, 62G07, 65D10, 68N15
h.ucv <- function(x,deriv.order=0,lower=0.1*hos,upper=2*hos,tol=0.1 * lower,
kernel="gaussian",...)
{
if (!is.numeric(x) || length(dim(x)) >=1 || length(x) < 3L)
stop("argument 'x' must be numeric and need at least 3 data points")
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
if (missing(kernel)) kernel <- "gaussian"
if (kernel != "gaussian") {stop("Expecting gaussian kernel.")}
name <- deparse(substitute(x))
x <- x[!is.na(x)]
x <- sort(x)
n <- length(x)
hos <- ((243 *(2*r+1)*A3_kMr(kernel,r))/(35* A2_kM(kernel)^2))^(1/(2*r+5)) * sd(x,na.rm = TRUE) * n^(-1/(2*r+5))
if (!is.numeric(upper)){
stop("argument 'upper' must be numeric. Default 2*hos (Oversmoothing) boundary was used")
upper= 2*hos
}
if (!is.numeric(lower)){
stop("argument 'lower' must be numeric. Default 0.1*hos boundary was used")
lower=0.1*hos
}
if (lower < 0 | lower >= upper){
stop("the boundaries must be positive and 'lower' must be smaller than 'upper'. Default boundaries were used")
upper=2*hos
lower=0.1*hos
}
R_Kr1 <- A3_kMr(kernel,r)
fucv <- function(h)
{
D <- kernel_fun_der(kernel, outer(x,x,"-")/h,deriv.order=2*r)
diag(D) <- 0
D <- ((-1)^r / ((n-1)*h^(2*r+1)))* colSums(D)
D1 <- mean(D)
D2 <- kernel_fun_conv(kernel,outer(x,x,"-")/h,deriv.order=r)
diag(D2) <- 0
D3 <- ((-1)^r / ((n-1)*h^(2*r+1)))* colSums(D2)
D4 <- mean(D3)
(1/(n*h^(2*r+1)))* R_Kr1 + D4 - 2*D1
}
obj <- optimize(fucv , c(lower, upper),tol=tol)
structure(list(x=x, data.name=name,n=n, kernel=kernel, deriv.order=r, h = obj$minimum ,
min.ucv=obj$objective),class="h.ucv")
}
A3_kMr <- function(kernel,r) {
xKr <- integrate(function(x) kernel_fun_der(kernel,x,deriv.order=r)^2, -Inf, Inf)$value
return(xKr)
}
A2_kM <-function(kernel) {
xKr <- 1
return(xKr)
}
####
#### r(th) derivative of Kernel functions K^r(x)
kernel_fun_der <- function(kernel,u,deriv.order=0)
{
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
Kr <- expression( dnorm(X) )
if (r == 0) {
DKr <- Kr
K <- function(X) eval(DKr);fx <- K(u)
} else {
if (r == 1){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)/sqrt(pi)}
else if (r == 2){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^2-1)/sqrt(pi) }
else if (r == 3){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^2-3)/sqrt(pi)}
else if (r == 4){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^4-6*u^2+3)/sqrt(pi)}
else if (r == 5){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^4-10*u^2+15)/sqrt(pi)}
else if (r == 6){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^6-15*u^4+45*u^2-15)/sqrt(pi)}
else if (r == 7){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^6-21*u^4+105*u^2-105)/sqrt(pi)}
else if (r == 8){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^8-28*u^6+210*u^4-420*u^2+105)/sqrt(pi)}
else if (r == 9){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^8-36*u^6+378*u^4-1260*u^2+945)/sqrt(pi)}
else if (r == 10){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^10-45*u^8+630*u^6-3150*u^4+4725*u^2-945)/sqrt(pi)}
else if (r == 11){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^10-55*u^8+990*u^6-6930*u^4+17325*u^2-10395)/sqrt(pi)}
else if (r == 12){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^12-66*u^10+1485*u^8-13860*u^6+51975*u^4-62370*u^2+10395)/sqrt(pi)}
else if (r == 13){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^12-78*u^10+2145*u^8-25740*u^6+135135*u^4-270270*u^2+135135)/sqrt(pi)}
else if (r == 14){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^14-91*u^12+3003*u^10-45045*u^8+315315*u^6-945945*u^4+945945*u^2-135135)/sqrt(pi)}
else if (r == 15){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^14-105*u^12+4095*u^10-75075*u^8+675675*u^6-2837835*u^4+4729725*u^2-2027025)/sqrt(pi)}
else if (r == 16){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^16-120*u^14+5460*u^12-120120*u^10+1351350*u^8-7567560*u^6+18918900*u^4-16216200*u^2+2027025)/sqrt(pi)}
else if (r == 17){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^16-136*u^14+7140*u^12-185640*u^10+2552550*u^8-18378360*u^6+64324260*u^4-91891800*u^2+34459425)/sqrt(pi)}
else if (r == 18){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^18-153*u^16+9180*u^14-278460*u^12+4594590*u^10-41351310*u^8+192972780*u^6-413513100*u^4+310134825*u^2-34459425)/sqrt(pi)}
else if (r == 19){fx <- -(1/2)*u*exp(-(1/2)*u^2)*sqrt(2)*(u^18-171*u^16+11628*u^14-406980*u^12+7936110*u^10-87297210*u^8+523783260*u^6-1571349780*u^4+1964187225*u^2-654729075)/sqrt(pi)}
else if (r == 20){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^20-190*u^18+14535*u^16-581400*u^14+13226850*u^12-174594420*u^10+1309458150*u^8-5237832600*u^6+9820936125*u^4-6547290750*u^2+654729075)/sqrt(pi)}
else if (r == 21){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^20-210*u^18+17955*u^16-813960*u^14+21366450*u^12-333316620*u^10+3055402350*u^8-15713497800*u^6+41247931725*u^4-45831035250*u^2+13749310575)/sqrt(pi)}
else if (r == 22){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^22-231*u^20+21945*u^18-1119195*u^16+33575850*u^14-611080470*u^12+6721885170*u^10-43212118950*u^8+151242416325*u^6-252070693875*u^4+151242416325*u^2-13749310575)/sqrt(pi)}
else if (r == 23){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^22-253*u^20+26565*u^18-1514205*u^16+51482970*u^14-1081142370*u^12+14054850810*u^10-110430970650*u^8+496939367925*u^6-1159525191825*u^4+1159525191825*u^2-316234143225)/sqrt(pi)}
else if (r == 24){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^24-276*u^22+31878*u^20-2018940*u^18+77224455*u^16-1853386920*u^14+28109701620*u^12-265034329560*u^10+1490818103775*u^8-4638100767300*u^6+6957151150950*u^4-3794809718700*u^2+316234143225)/sqrt(pi)}
else if (r == 25){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^24-300*u^22+37950*u^20-2656500*u^18+113565375*u^16-3088978200*u^14+54057118500*u^12-602350749000*u^10+4141161399375*u^8-16564645597500*u^6+34785755754750*u^4-31623414322500*u^2+7905853580625)/sqrt(pi)}
else if (r == 26){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^26-325*u^24+44850*u^22-3453450*u^20+164038875*u^18-5019589575*u^16+100391791500*u^14-1305093289500*u^12+10767019638375*u^10-53835098191875*u^8+150738274937250*u^6-205552193096250*u^4+102776096548125*u^2-7905853580625)/sqrt(pi)}
else if (r == 27){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^26-351*u^24+52650*u^22-4440150*u^20+233107875*u^18-7972289325*u^16+180705224700*u^14-2710578370500*u^12+26428139112375*u^10-161505294575625*u^8+581419060472250*u^6-1109981842719750*u^4+924984868933125*u^2-213458046676875)/sqrt(pi)}
else if (r == 28){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^28-378*u^26+61425*u^24-5651100*u^22+326351025*u^20-12401338950*u^18+316234143225*u^16-5421156741000*u^14+61665657928875*u^12-452214824811750*u^10+2034966711652875*u^8-5179915266025500*u^6+6474894082531875*u^4-2988412653476250*u^2+213458046676875)/sqrt(pi)}
else if (r == 29){fx <- -(1/2)*exp(-(1/2)*u^2)*sqrt(2)*u*(u^28-406*u^26+71253*u^24-7125300*u^22+450675225*u^20-18928359450*u^18+539458244325*u^16-10480903032600*u^14+137561852302875*u^12-1192202719958250*u^10+6557114959770375*u^8-21459648959248500*u^6+37554385678684875*u^4-28887988983603750*u^2+6190283353629375)/sqrt(pi)}
else if (r == 30){fx <- (1/2)*exp(-(1/2)*u^2)*sqrt(2)*(u^30-435*u^28+82215*u^26-8906625*u^24+614557125*u^22-28392539175*u^20+899097073875*u^18-19651693186125*u^16+294775397791875*u^14-2980506799895625*u^12+19671344879311125*u^10-80473683597181875*u^8+187771928393424375*u^6-216659917377028125*u^4+92854250304440625*u^2-6190283353629375)/sqrt(pi)}
else {K <- function(X) (-1)^r * Hermite(X,r) * eval(Kr);fx <- K(u)}
}
return(fx)
}
#### Hermite Polynomial
Hermite <-function (x, n, prob = TRUE) {
if (any(n < 0 || n != round(n))) stop("Argument 'n' must be a vector of non-negative integers")
if ((length(n) != 1) && (length(x) != length(n)) && (length(x) != 1))
stop(paste("Argument 'n' must be either a vector of same length",
"as argument 'x',\n a single integer or 'x' must be a ",
"single value!", sep = ""))
H <- function(x, n) {
if (n <= 1) {
return(switch(n + 1, 1, x))
}
else {
return(x * Recall(x, n - 1) - (n - 1) * Recall(x, n - 2))
}
}
scale <- 1
if (!prob) {
x <- sqrt(2) * x
scale <- 2^(n/2)
}
scale * mapply(H, x, n)
}
#### Kernels density convolution
kernel_fun_conv <- function(kernel,u,deriv.order=0)
{
if (any(deriv.order < 0 || deriv.order != round(deriv.order)))
stop("argument 'deriv.order' is non-negative integers")
r <- deriv.order
if (r==0) {fx <- dnorm(u,mean=0,sd=sqrt(2))}
else if (r==1) {fx <- (1/8)*exp(-(1/4)*u^2)*(u^2-2)/sqrt(pi)}
else if (r==2) {fx <- (1/32)*exp(-(1/4)*u^2)*(12-12*u^2+u^4)/sqrt(pi)}
else if (r==3) {fx <- (1/128)*exp(-(1/4)*u^2)*(u^6-30*u^4+180*u^2-120)/sqrt(pi)}
else if (r==4) {fx <- (1/512)*exp(-(1/4)*u^2)*(u^8-56*u^6+840*u^4-3360*u^2+1680)/sqrt(pi)}
else if (r==5) {fx <- (1/2048)*exp(-(1/4)*u^2)*(u^10-90*u^8+2520*u^6-25200*u^4+75600*u^2-30240)/sqrt(pi)}
else if (r==6) {fx <- (1/8192)*exp(-(1/4)*u^2)*(u^12-132*u^10+5940*u^8-110880*u^6+831600*u^4-1995840*u^2+665280)/sqrt(pi)}
else if (r==7) {fx <- (1/32768)*exp(-(1/4)*u^2)*(u^14-182*u^12+12012*u^10-360360*u^8+5045040*u^6-30270240*u^4+60540480*u^2-17297280)/sqrt(pi)}
else if (r==8) {fx <- (1/131072)*exp(-(1/4)*u^2)*(u^16-240*u^14+21840*u^12-960960*u^10+21621600*u^8-242161920*u^6+1210809600*u^4-2075673600*u^2+518918400)/sqrt(pi)}
else if (r==9) {fx <- (1/524288)*exp(-(1/4)*u^2)*(u^18-306*u^16+36720*u^14-2227680*u^12+73513440*u^10-1323241920*u^8+12350257920*u^6-52929676800*u^4+79394515200*u^2-17643225600)/sqrt(pi)}
else if (r==10) {fx <- (1/2097152)*exp(-(1/4)*u^2)*(u^20-380*u^18+58140*u^16-4651200*u^14+211629600*u^12-5587021440*u^10+83805321600*u^8-670442572800*u^6+2514159648000*u^4-3352212864000*u^2+670442572800)/sqrt(pi)}
else if (r==11) {fx <- (1/8388608)*exp(-(1/4)*u^2)*(u^22-462*u^20+87780*u^18-8953560*u^16+537213600*u^14-19554575040*u^12+430200650880*u^10-5531151225600*u^8+38718058579200*u^6-129060195264000*u^4+154872234316800*u^2-28158588057600)/sqrt(pi)}
else if (r==12) {fx <- (1/33554432)*exp(-(1/4)*u^2)*(u^24-552*u^22+127512*u^20-16151520*u^18+1235591280*u^16-59308381440*u^14+1799020903680*u^12-33924394183680*u^10+381649434566400*u^8-2374707592857600*u^6+7124122778572800*u^4-7771770303897600*u^2+1295295050649600)/sqrt(pi)}
else if (r==13) {fx <- (1/134217728)*exp(-(1/4)*u^2)*(u^26-650*u^24+179400*u^22-27627600*u^20+2624622000*u^18-160626866400*u^16+6425074656000*u^14-167051941056000*u^12+2756357027424000*u^10-27563570274240000*u^8+154355993535744000*u^6-420970891461120000*u^4+420970891461120000*u^2-64764752532480000)/sqrt(pi)}
else if (r==14) {fx <- (1/536870912)*exp(-(1/4)*u^2)*(u^28-756*u^26+245700*u^24-45208800*u^22+5221616400*u^20-396842846400*u^18+20238985166400*u^16-693908062848000*u^14+15786408429792000*u^12-231533990303616000*u^10+2083805912732544000*u^8-10608466464820224000*u^6+26521166162050560000*u^4-24481076457277440000*u^2+3497296636753920000)/sqrt(pi)}
else if (r==15) {fx <- (1/2147483648)*exp(-(1/4)*u^2)*(u^30-870*u^28+328860*u^26-71253000*u^24+9832914000*u^22-908561253600*u^20+57542212728000*u^18-2515416727824000*u^16+75462501834720000*u^14-1526019481546560000*u^12+20143457156414592000*u^10-164810104007028480000*u^8+769113818699466240000*u^6-1774878043152614400000*u^4+1521324036987955200000*u^2-202843204931727360000)/sqrt(pi)}
else if (r==16) {fx <- (1/8589934592)*exp(-(1/4)*u^2)*(u^32-992*u^30+431520*u^28-108743040*u^26+17670744000*u^24-1950850137600*u^22+150215460595200*u^20-8154553575168000*u^18+311911674250176000*u^16-8317644646671360000*u^14+151381132569418752000*u^12-1816573590833025024000*u^10+13624301931247687680000*u^8-58689300626913116160000*u^6+125762787057670963200000*u^4-100610229646136770560000*u^2+12576278705767096320000)/sqrt(pi)}
else if (r==17) {fx <- (1/34359738368)*exp(-(1/4)*u^2)*(u^34-1122*u^32+556512*u^30-161388480*u^28+30502422720*u^26-3965314953600*u^24+364808975731200*u^22-24077392398259200*u^20+1143676138917312000*u^18-38884988723188608000*u^16+933239729356526592000*u^14-15440875522080712704000*u^12+169849630742887839744000*u^10-1175882058989223505920000*u^8+4703528235956894023680000*u^6-9407056471913788047360000*u^4+7055292353935341035520000*u^2-830034394580628357120000)/sqrt(pi)}
else if (r==18) {fx <- (1/137438953472)*exp(-(1/4)*u^2)*(u^36-1260*u^34+706860*u^32-233735040*u^30+50837371200*u^28-7686610525440*u^26+832716140256000*u^24-65665615631616000*u^22+3792189302725824000*u^20-160114659448423680000*u^18+4899508579121764608000*u^16-106898368999020318720000*u^14+1621291929818474833920000*u^12-16462348825849129082880000*u^10+105829385309030115532800000*u^8-395096371820379097989120000*u^6+740805697163210808729600000*u^4-522921668585795864985600000*u^2+58102407620643984998400000)/sqrt(pi)}
else if (r==19) {fx <- (1/549755813888)*exp(-(1/4)*u^2)*(u^38-1406*u^36+885780*u^34-331281720*u^32+82157866560*u^30-14295468781440*u^28+1801229066461440*u^26-167256984742848000*u^24+11540731947256512000*u^22-592424239959167616000*u^20+22512121118448369408000*u^18-626246278385927367168000*u^16+12524925567718547343360000*u^14-175348957948059662807040000*u^12+1653290174938848249323520000*u^10-9919741049633089495941120000*u^8+34719093673715813235793920000*u^6-61268988835969082180812800000*u^4+40845992557312721453875200000*u^2-4299578163927654889881600000)/sqrt(pi)}
else if (r==20) {fx <- (1/2199023255552)*exp(-(1/4)*u^2)*(u^40-1560*u^38+1096680*u^36-460605600*u^34+129199870800*u^32-25633254366720*u^30+3716821883174400*u^28-401416763382835200*u^26+32615112024855360000*u^24-2000393537524462080000*u^22+92418181433630148096000*u^20-3192628085889041479680000*u^18+81412016190170557731840000*u^16-1502991068126225681203200000*u^14+19538883885640933855641600000*u^12-171942178193640217929646080000*u^10+967174752339226225854259200000*u^8-3185987419470392273402265600000*u^6+5309979032450653789003776000000*u^4-3353670967863570814107648000000*u^2+335367096786357081410764800000)/sqrt(pi)}
else if (r==21) {fx <- (1/8796093022208)*exp(-(1/4)*u^2)*(u^42-1722*u^40+1343160*u^38-629494320*u^36+198290710800*u^34-44496435503520*u^32+7356744003248640*u^30-914338183260902400*u^28+86404958318155276800*u^26-6240358100755658880000*u^24+344467767161712370176000*u^22-14467646220791919547392000*u^20+458142130325077452334080000*u^18-10783960913805669262632960000*u^16+184867901379525758787993600000*u^14-2243063870071579206627655680000*u^12+18505276928090528454678159360000*u^10-97969113148714562407119667200000*u^8+304792796462667527488816742400000*u^6-481251783888422411824447488000000*u^4+288751070333053447094668492800000*u^2-27500101936481280675682713600000)/sqrt(pi)}
else if (r==22) {fx <- (1/35184372088832)*exp(-(1/4)*u^2)*(u^44-1892*u^42+1629012*u^40-847086240*u^38+297750813360*u^36-75033204966720*u^34+14031209328776640*u^32-1988422807735203840*u^30+216240980341203417600*u^28-18164242348661087078400*u^26+1180675752662970660096000*u^24-59248455951814527670272000*u^22+2281065554144859315305472000*u^20-66677300813465118447390720000*u^18+1457375289208594731778682880000*u^16-23318004627337515708458926080000*u^14+265242302635964241183720284160000*u^12-2059528467526310578603004559360000*u^10+10297642337631552893015022796800000*u^8-30350945837229840105728488243200000*u^6+45526418755844760158592732364800000*u^4-26015096431911291519195847065600000*u^2+2365008766537390138108713369600000)/sqrt(pi)}
else if (r==23) {fx <- (1/140737488355328)*exp(-(1/4)*u^2)*(u^46-2070*u^44+1958220*u^42-1124018280*u^40+438367129200*u^38-123268836731040*u^36+25886455713518400*u^34-4149229044366806400*u^32+514504401501483993600*u^30-49735425478476786048000*u^28+3759998166172845025228800*u^26-222181709819304478763520000*u^24+10220358651688006023121920000*u^22-363215822852296829437102080000*u^20+9858715191705199656149913600000*u^18-201117789910786072985458237440000*u^16+3016766848661791094781873561600000*u^14-32297150968026234073547116953600000*u^12+236845773765525716539345524326400000*u^10-1121901033626174446765320904704000000*u^8+3141322894153288450942898533171200000*u^6-4487604134504697787061283618816000000*u^4+2447784073366198792942518337536000000*u^2-212850788988365112429784203264000000)/sqrt(pi)}
else if (r==24) {fx <- (1/562949953421312)*exp(-(1/4)*u^2)*(u^48-2256*u^46+2334960*u^44-1472581440*u^42+633946309920*u^40-197791248695040*u^38+46349082610871040*u^36-8342834869956787200*u^34+1170082590511439404800*u^32-128969103309705321062400*u^30+11220311987944362932428800*u^28-771141442080539852446924800*u^26+41770161446029242007541760000*u^24-1773625316785241660627927040000*u^22+58529635453912974800721592320000*u^20-1482750764832462028284947005440000*u^18+28357608377420836290949611479040000*u^16-400342706504764747636935691468800000*u^14+4047909587992621337217905324851200000*u^12-28122319242896106132250710677913600000*u^10+126550436593032477595128198050611200000*u^8-337467830914753273587008528134963200000*u^6+460183405792845373073193447456768000000*u^4-240095689978875846820796581281792000000*u^2+20007974164906320568399715106816000000)/sqrt(pi)}
else if (r==25) {fx <- (1/2251799813685248)*exp(-(1/4)*u^2)*(u^50-2450*u^48+2763600*u^46-1906884000*u^44+901956132000*u^42-310633691860800*u^40+80764759883808000*u^38-16222178913804864000*u^36+2554993178924266080000*u^34-318522482972558504640000*u^32+31597430310877803660288000*u^30-2499069488223971744040960000*u^28+157441377758110219874580480000*u^26-7872068887905510993729024000000*u^24+310384430437417290609887232000000*u^22-9559840457472452550784526745600000*u^20+227046210864970748081132510208000000*u^18-4086831795569473465460385183744000000*u^16+54491090607592979539471802449920000000*u^14-521967288977995909272835160309760000000*u^12+3444984107254773001200712058044416000000*u^10-14764217602520455719431623105904640000000*u^8+37581644806415705467644131542302720000000*u^6-49019536704020485392579302011699200000000*u^4+24509768352010242696289651005849600000000*u^2-1960781468160819415703172080467968000000)/sqrt(pi)}
else if (r==26) {fx <- (1/9007199254740992)*exp(-(1/4)*u^2)*(u^52-2652*u^50+3248700*u^48-2443022400*u^46+1264264092000*u^44-478397532412800*u^42+137300091802473600*u^40-30598306173122688000*u^38+5377652309926312416000*u^36-752871323389683738240000*u^34+84472162484322515430528000*u^32-7617853198586175937007616000*u^30+552294356897497755433052160000*u^28-32118041062654484854414417920000*u^26+1491194763623243939669240832000000*u^24-54875967301335376979828062617600000*u^22+1584543555826059010292535308083200000*u^20-35419208894935436700656671592448000000*u^18+602126551213902423911163417071616000000*u^16-7605809067965083249404169478799360000000*u^14+69212862518482257569577942257074176000000*u^12-435052278687602761865918494187323392000000*u^10+1779759321903829480360575658039050240000000*u^8-4333327044635410908704010297834209280000000*u^6+5416658805794263635880012872292761600000000*u^4-2599996226781246545222406178700525568000000*u^2+199999709752403580401723552207732736000000)/sqrt(pi)}
else if (r==27) {fx <- (1/36028797018963968)*exp(-(1/4)*u^2)*(u^54-2862*u^52+3795012*u^50-3099259800*u^48+1747982527200*u^46-723664766260800*u^44+228195622960905600*u^42-56136123248382777600*u^40+10946544033434641632000*u^38-1710093434556567348288000*u^36+215471772754127485884288000*u^34-21978120820921003560197376000*u^32+1816857987862802960976316416000*u^30-121589726880049121234568867840000*u^28+6565845251522652546666718863360000*u^26-284519960899314943688891150745600000*u^24+9815938651026365557266744700723200000*u^22-266762568045540052203366826572595200000*u^20+5631654214294734435404410783199232000000*u^18-90699273135483617749144721034682368000000*u^16+1088391277625803412989736652416188416000000*u^14-9432724406090296245911050987606966272000000*u^12+56596346436541777475466305925641797632000000*u^10-221463964316902607512694240578598338560000000*u^8+516749250072772750862953228016729456640000000*u^6-620099100087327301035543873620075347968000000*u^4+286199584655689523554866403209265545216000000*u^2-21199969233754779522582696534019670016000000)/sqrt(pi)}
else if (r==28) {fx <- (1/144115188075855872)*exp(-(1/4)*u^2)*(u^56-3080*u^54+4407480*u^52-3896212320*u^50+2386430046000*u^48-1076757236755200*u^46+371481246680544000*u^44-100406074102798464000*u^42+21612407450627369376000*u^40-3746150624775410691840000*u^38+526708777843422743272704000*u^36-60332096371155696047600640000*u^34+5641051010703057580450659840000*u^32-430455584816725624600542658560000*u^30+26749739913610806671605150924800000*u^28-1348186891645984656248899606609920000*u^26+54770092473118126660111546518528000000*u^24-1778417120303600348022445510483968000000*u^22+45646039421125742265909434769088512000000*u^20-912920788422514845318188695381770240000000*u^18+13967688062864477133368287039341084672000000*u^16-159630720718451167238494709021040967680000000*u^14+1320581416852641474427547138264975278080000000*u^12-7578989001067333679323314010912032030720000000*u^10+28421208754002501297462427540920120115200000000*u^8-63663507608965602906315837691661069058048000000*u^6+73457893394960311045749043490378156605440000000*u^4-32647952619982360464777352662390291824640000000*u^2+2331996615713025747484096618742163701760000000)/sqrt(pi)}
else if (r==29) {fx <- (1/576460752303423488)*exp(-(1/4)*u^2)*(u^58-3306*u^56+5091240*u^54-4857042960*u^52+3220219482480*u^50-1577907546415200*u^48+593293237452115200*u^46-175445285932268352000*u^44+41492810122981465248000*u^42-7938957670197120350784000*u^40+1238477396550750774722304000*u^38-158299929050032326296323584000*u^36+16621492550253394261113976320000*u^34-1434562664721869873920760110080000*u^32+101649011671721065352099573514240000*u^30-5895642676959821790421775263825920000*u^28+278569116486351579597428881215774720000*u^26-10651172100948736866960516046485504000000*u^24+326635944429094597253455825425555456000000*u^22-7942410859275879154268241649821401088000000*u^20+150905806326241703931096591346606620672000000*u^18-2198913177896664828710264616764839329792000000*u^16+23988143758872707222293795819252792688640000000*u^14-189819224526731857150324819091478620405760000000*u^12+1044005734897025214326786505003132412231680000000*u^10-3758420645629290771576431418011276684034048000000*u^8+8095059852124626277241544592639672857919488000000*u^6-8994510946805140308046160658488525397688320000000*u^4+3854790405773631560591211710780796599009280000000*u^2-265847614191284935213187014536606662000640000000)/sqrt(pi)}
else if (r==30) {fx <- (1/2305843009213693952)*exp(-(1/4)*u^2)*(u^60-3540*u^58+5851620*u^56-6007663200*u^54+4298483019600*u^52-2279915393595840*u^50+930965452384968000*u^48-300036865797212544000*u^46+77634539025028745760000*u^44-16320505315039376330880000*u^42+2810391015249780604177536000*u^40-398564543980877976592450560000*u^38+46698479069759536257415457280000*u^36-4526160279069001206487959705600000*u^34+362739416651101382405677913548800000*u^32-23989166754526171423095499349360640000*u^30+1304410942277360571130817777121484800000*u^28-58007921903628505398523425853167206400000*u^26+2094730513186584917168901489142149120000000*u^24-60857433856789203909328085368761384960000000*u^22+1405806722091830610305478772018387992576000000*u^20-25438407352137887234099139684142258913280000000*u^18+353825120443372431528833488333978692157440000000*u^16-3692088213322147111605219008702386352947200000000*u^14+27998335617692948929672910815993096509849600000000*u^12-147831212061418770348672969108443549572005888000000*u^10+511723426366449589668483354606150748518481920000000*u^8-1061352291723006556349446957701645996927221760000000*u^6+1137163169703221310374407454680334996707737600000000*u^4-470550277118574335327341015729793791741132800000000*u^2+31370018474571622355156067715319586116075520000000)/sqrt(pi)}
else if (r>=31) {fx <- NA}
return(fx)
}
Try the shazam package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.