R/utils.R

## Sat Jul 20 02:11:02 2013
## Original file Copyright 2013 A.C. Guidoum
## This file is part of the R package kedd.
## Arsalane Chouaib GUIDOUM <acguidoum@usthb.dz> and <starsalane@gmail.com> 
## Department of Probabilities-Statistics
## Faculty of Mathematics
## University of Science and Technology Houari Boumediene
## BP 32 El-Alia, U.S.T.H.B, Algeris
## Algeria

## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## A copy of the GNU General Public License is available at
## http://www.r-project.org/Licenses/
## Unlimited use and distribution (see LICENCE).
## kedd : Kernel estimator and bandwidth selection for density and its derivatives.
###################################################################################################



####
#### 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
     if (kernel=="gaussian")          {Kr <- expression( dnorm(X) )}
     else if (kernel=="epanechnikov") {Kr <- expression( (3/4)*(1-X^2)) }
     else if (kernel=="uniform")      {Kr <- expression( 0.5 ) }
     else if (kernel=="triweight")    {Kr <- expression( (35/32)*(1-X^2)^3 )}
     else if (kernel=="biweight")     {Kr <- expression( (15/16)*(1-X^2)^2 )}
     else if (kernel=="cosine")       {Kr <- expression( (pi/4)*cos((pi*X)/2) )}
     else if (kernel=="triangular")   {Kr <- expression(1-abs(X));Kr1 <- expression(1-X); Kr2 <- expression(1+X)}
     else if (kernel=="tricube")      {Kr <- expression((70/81)*(1-abs(X)^3)^3);Kr1 <- expression((70/81)*(1-X^3)^3); Kr2 <- expression((70/81)*(1+X^3)^3)}
     else if (kernel=="silverman")    {r <- r%%8;Kr <- expression(0.5*exp(-abs(X)/sqrt(2))*sin((abs(X)/sqrt(2))+0.25*pi));Kr1 <- expression(0.5*exp(-X/sqrt(2)) * sin((X/sqrt(2)) + 0.25*pi)); Kr2 <- expression(0.5*exp(X/sqrt(2)) * sin((-X/sqrt(2)) + 0.25*pi))}

     if (kernel=="epanechnikov" && r >= 3) stop(" 'epanechnikov kernel derivative = 0' for 'order >= 3' ")
     if (kernel=="uniform" && r >= 1)      stop(" 'uniform kernel derivative = 0' for 'order >= 1' ")
     if (kernel=="triweight" && r >= 7)    stop(" 'triweight kernel derivative = 0' for 'order >= 7' ")
     if (kernel=="biweight" && r >= 5)     stop(" 'biweight kernel derivative = 0' for 'order >= 5' ")
     if (kernel=="triangular" && r >= 2)   stop(" 'triangular kernel derivative = 0' for 'order >= 2' ")
     if (kernel=="tricube" && r >= 10)     stop(" 'tricube kernel derivative = 0' for 'order >= 10' ")

        if (r == 0) {DKr <- Kr
        if (kernel=="gaussian" || kernel =="silverman"){
             K <- function(X) eval(DKr);fx <- K(u)}else{
             K <- function(X) eval(DKr)* (X >= -1 & X <= 1);fx <- K(u)}
        } else { if (kernel=="gaussian"){
		     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)}}	
        else if (kernel=="cosine"){
		     if (r%%2==0){
			 if ((r%/%2)%%2==0){
			 K <- function(X)( ((-1)^r/2^(r+2))*pi^(r+1)*cos(0.5*pi*X))* (X >= -1 & X <= 1)}else{
			 K <- function(X)(-((-1)^r/2^(r+2))*pi^(r+1)*cos(0.5*pi*X))* (X >= -1 & X <= 1)}}else{
			 if ((r%/%2)%%2==0){
			 K <- function(X)( ((-1)^r/2^(r+2))*pi^(r+1)*sin(0.5*pi*X))* (X >= -1 & X <= 1)}else{
			 K <- function(X)(-((-1)^r/2^(r+2))*pi^(r+1)*sin(0.5*pi*X))* (X >= -1 & X <= 1)}}
             fx <- K(u)}
        else if (kernel=="tricube"){
             if (r == 1){fx <- (u < 0 & u >= -1) *((70/9)*(u^3+1)^2*u^2)+ (u >= 0 & u <= 1)*(-(70/9)*(u^3-1)^2*u^2)}
		     else if (r == 2){fx <- (u < 0 & u >= -1) *((560/9)*u^7+(700/9)*u^4+(140/9)*u)+ (u >= 0 & u <= 1)*(-(560/9)*u^7+(700/9)*u^4-(140/9)*u)}
		     else if (r == 3){fx <- (u < 0 & u >= -1) *((3920/9)*u^6+(2800/9)*u^3+(140/9))+ (u >= 0 & u <= 1)*(-(3920/9)*u^6+(2800/9)*u^3-(140/9))}
		     else if (r == 4){fx <- (u < 0 & u >= -1) *((7840/3)*u^5+(2800/3)*u^2)+ (u >= 0 & u <= 1)*(-(7840/3)*u^5+(2800/3)*u^2)}
		     else if (r == 5){fx <- (u < 0 & u >= -1) *((39200/3)*u^4+(5600/3)*u)+ (u >= 0 & u <= 1)*(-(39200/3)*u^4+(5600/3)*u)}
		     else if (r == 6){fx <- (u < 0 & u >= -1) *((156800/3)*u^3+(5600/3))+ (u >= 0 & u <= 1)*(-(156800/3)*u^3+(5600/3))}
		     else if (r == 7){fx <- (u <0 & u >= -1) *(156800*u^2)+ (u >= 0 & u <= 1)*(-156800*u^2)}
		     else if (r == 8){fx <- (u < 0 & u >= -1) *(313600*u)+ (u >= 0 & u <= 1)*( -313600*u)}
		     else if (r == 9){fx <- (u < 0 & u >= -1) *(313600)+ (u >= 0 & u <= 1)*( -313600)}}
        else if (kernel=="triweight"){
             if (r == 1){fx <- (-(105/16)*(u^2-1)^2*u)*(u >= -1 & u <= 1)}
		     else if (r == 2){fx <- (-(525/16)*u^4+(315/8)*u^2-(105/16))*(u >= -1 & u <= 1)}
		     else if (r == 3){fx <- (-(525/4)*u^3+(315/4)*u)*(u >= -1 & u <= 1)}
		     else if (r == 4){fx <- (-(1575/4)*u^2+(315/4))*(u >= -1 & u <= 1)}
		     else if (r == 5){fx <- (-(1575/2)*u)*(u >= -1 & u <= 1)}
		     else if (r == 6){fx <- (-1575/2)*(u >= -1 & u <= 1)}}
        else if (kernel=="biweight"){
             if (r == 1){fx <- (((15/4)*(u^2-1))*u)*(u >= -1 & u <= 1)}
		     else if (r == 2){fx <- ((45/4)*u^2-(15/4))*(u >= -1 & u <= 1)}
		     else if (r == 3){fx <- ((45/2)*u)*(u >= -1 & u <= 1)}
		     else if (r == 4){fx <- (45/2)*(u >= -1 & u <= 1)}}	
        else if (kernel =="triangular"){
             if (r == 1){fx <- (u <= 0 & u >= -1) - (u >= 0 & u <= 1)}}
        else if (kernel=="epanechnikov"){
             if (r == 1){fx <- (-(3/2)*u)*(u >= -1 & u <= 1)}
		     else if (r == 2){fx <- (-3/2)*(u >= -1 & u <= 1)}}				 
        else if (kernel=="silverman"){
             if (r == 1){fx <- (u < 0)*((1/4)*sqrt(2)*exp((1/2)*u*sqrt(2))*(sin(-(1/2)*u*sqrt(2)+(1/4)*pi)-cos(-(1/2)*u*sqrt(2)+(1/4)*pi)))+ (u >= 0)* (-(1/4)*sqrt(2)*exp(-(1/2)*u*sqrt(2))*(sin((1/2)*u*sqrt(2)+(1/4)*pi)-cos((1/2)*u*sqrt(2)+(1/4)*pi)))}
		     else if (r == 2){fx <- (u < 0)*(-(1/2)*exp((1/2)*u*sqrt(2))*cos(-(1/2)*u*sqrt(2)+(1/4)*pi))+ (u >= 0)* (-(1/2)*exp(-(1/2)*u*sqrt(2))*cos((1/2)*u*sqrt(2)+(1/4)*pi))}
		     else if (r == 3){fx <- (u < 0)*(-(1/4)*exp((1/2)*u*sqrt(2))*sqrt(2)*(cos(-(1/2)*u*sqrt(2)+(1/4)*pi)+sin(-(1/2)*u*sqrt(2)+(1/4)*pi)))+ (u >= 0)* ((1/4)*exp(-(1/2)*u*sqrt(2))*sqrt(2)*(cos((1/2)*u*sqrt(2)+(1/4)*pi)+sin((1/2)*u*sqrt(2)+(1/4)*pi)))}
		     else if (r == 4){fx <- (u < 0)*(-(1/2)*exp((1/2)*u*sqrt(2))*sin(-(1/2)*u*sqrt(2)+(1/4)*pi))+ (u >= 0)* (-(1/2)*exp(-(1/2)*u*sqrt(2))*sin((1/2)*u*sqrt(2)+(1/4)*pi))}
		     else if (r == 5){fx <- (u < 0)*(-(1/4)*sqrt(2)*exp((1/2)*u*sqrt(2))*(sin(-(1/2)*u*sqrt(2)+(1/4)*pi)-cos(-(1/2)*u*sqrt(2)+(1/4)*pi)))+ (u >= 0)* ((1/4)*sqrt(2)*exp(-(1/2)*u*sqrt(2))*(sin((1/2)*u*sqrt(2)+(1/4)*pi)-cos((1/2)*u*sqrt(2)+(1/4)*pi)))}
		     else if (r == 6){fx <- (u < 0)*((1/2)*exp((1/2)*u*sqrt(2))*cos(-(1/2)*u*sqrt(2)+(1/4)*pi))+ (u >= 0)* ((1/2)*exp(-(1/2)*u*sqrt(2))*cos((1/2)*u*sqrt(2)+(1/4)*pi))}
		     else if (r == 7){fx <- (u < 0)*((1/4)*exp((1/2)*u*sqrt(2))*sqrt(2)*(cos(-(1/2)*u*sqrt(2)+(1/4)*pi)+sin(-(1/2)*u*sqrt(2)+(1/4)*pi)))+ (u >= 0)* (-(1/4)*exp(-(1/2)*u*sqrt(2))*sqrt(2)*(cos((1/2)*u*sqrt(2)+(1/4)*pi)+sin((1/2)*u*sqrt(2)+(1/4)*pi)))}}
               }
return(fx)
}


####
####  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 (kernel=="epanechnikov" && r >= 3) stop(" 'epanechnikov kernel derivative = 0' for 'order >= 3' ")
   if (kernel=="uniform" && r >= 1)      stop(" 'uniform kernel derivative = 0' for 'order >= 1' ")
   if (kernel=="triweight" && r >= 7)    stop(" 'triweight kernel derivative = 0' for 'order >= 7' ")
   if (kernel=="biweight" && r >= 5)     stop(" 'biweight kernel derivative = 0' for 'order >= 5' ")
   if (kernel=="triangular" && r >= 2)   stop(" 'triangular kernel derivative = 0' for 'order >= 2' ")
   if (kernel=="tricube" && r >= 10)     stop(" 'tricube kernel derivative = 0' for 'order >= 10' ")

   if (kernel=="uniform"){
       if (r==0){fx <- (0.5-0.25*abs(u))*(abs(u)<=2)}
	}
   else if(kernel=="epanechnikov"){
       if (r==0){fx <- ((3/5)-(3/160)*abs(u)^5+(3/8)*abs(u)^3-(3/4)*u^2)*(abs(u)<=2)}
	   else if (r==1){fx <- (-(3/2)-(3/8)*abs(u)^3+(9/4)*abs(u))*(abs(u)<=2)}
	   else if (r==2){fx <- ((9/2)-(9/4)*abs(u))*(abs(u)<=2)}
   }
   else if (kernel=="triweight"){
       if (r==0) {fx <- (-(35/22)*u^2+(350/429)+(35/11264)*abs(u)^11-(35/768)*abs(u)^9+
	                     (35/64)*abs(u)^7-(35/32)*u^6+(35/24)*u^4-(175/1757184)*abs(u)^13)*
						 (abs(u)<=2)}
	   else if (r==1) {fx <- (-(35/11)+(735/32)*abs(u)^5+(35/2)*u^2-(175/11264)*abs(u)^11+
	                     (175/512)*abs(u)^9-(105/32)*abs(u)^7-(525/16)*u^4)*(abs(u)<=2)}
	   else if (r==2) {fx <- (35-(2205/16)*abs(u)^5+(3675/8)*abs(u)^3-(1575/4)*u^2-(875/512)*
	                     abs(u)^9+(1575/64)*abs(u)^7)*(abs(u)<=2)}
	   else if (r==3) {fx <- ((11025/4)*abs(u)+(33075/32)*abs(u)^5-(11025/4)*abs(u)^3-(7875/64)*
	                     abs(u)^7-(1575/2))*(abs(u)<=2)}
	   else if (r==4) {fx <- (33075-(165375/32)*abs(u)^5+(165375/8)*abs(u)^3+(165375/4)*u^2-
	                     99225*abs(u))*(abs(u)<=2)}
	   else if (r==5) {fx <- (-(826875/2)-(826875/8)*abs(u)^3+(2480625/4)*abs(u))*(abs(u)<=2)}
	   else if (r==6) {fx <- ((2480625/2)-(2480625/4)*abs(u))*(abs(u)<=2)}
   }
   else if (kernel=="biweight"){
       if (r==0) {fx <- (-(15/32)*abs(u)^5-(15/14)*u^2-(5/3584)*abs(u)^9+(15/448)*abs(u)^7+
	                    (15/16)*u^4+(5/7))*(abs(u)<=2)}
	   else if (r==1) {fx <- ((45/32)*abs(u)^5-(75/8)*abs(u)^3+(45/4)*u^2-(45/448)*abs(u)^7-
	                     (15/7)) *(abs(u)<=2)}
	   else if (r==2) {fx <- ((45/2)-(135/32)*abs(u)^5+(225/8)*abs(u)^3-(225/4)*abs(u)) *
	                     (abs(u)<=2)}
	   else if (r==3) {fx <- (-(675/2)-(675/8)*abs(u)^3+(2025/4)*abs(u)) *(abs(u)<=2)}
	   else if (r==4) {fx <- ((2025/2)-(2025/4)*abs(u)) *(abs(u)<=2)}
   }
   else if (kernel=="cosine"){
       fx <- ((-1)^(r+1) * (pi^(2*r+2)/2^(2*r+5)) * cos(0.5*pi*abs(u))*abs(u) + (-1)^r * (pi^(2*r+2)/
	   2^(2*r+4))* cos(0.5*pi*abs(u))+(pi^(2*r+1)/2^(2*r+4))*sin(0.5*pi*abs(u)))*(abs(u)<=2)
   }
   else if (kernel=="triangular"){
       if (r==0) {fx <- ((1/6)*(u+2)^3 *(u <=- 1 & u > -2 ) + (-(1/2)*u^3-u^2+(2/3)) *
	                  (u <= 0 & u > -1 )+((2/3)+(1/2)*u^3-u^2)*(u <= 1 & u > 0 )-(1/6)*(u-2)^3 *
					  (u <=2 & u >1) )}
	   else if (r==1) {fx <- (u+2)*(u <=-1 & u >= -2 )+(-3*u-2)*(u <=0 & u > -1 )+
                         (3*u-2)*(u <= 1 & u > 0 )+(2-u)*(u <= 2 & u > 1)}   
   }
   else if (kernel=="tricube"){
       if (r==0) {fx <- (u <= -1 & u > -2)*((2870/351)*u^6+(21560/2187)*u^4+(7840/1683)*u^2+
	                  (22400/20007)+(70/81)*u^9+(15680/6561)*u+(11305/2187)*u^7+(980/99)*u^5+
					  (2800/351)*u^3+(35/625482)*u^16+(1085/6561)*u^10+(245/303046029)*u^19+
					  (665/312741)*u^13+(245/99)*u^8)+ (u <= 0 & u > -1)*(-(350/117)*u^6+(980/729)*u^4-
					  (210/187)*u^2+(175/247)-(70/81)*u^9-(2905/729)*u^7-(35/625482)*u^16-
					  (1085/6561)*u^10-(245/101015343)*u^19-(1295/312741)*u^13-(245/99)*u^8)+
					  (u > 0 & u <= 1)*(-(350/117)*u^6+(980/729)*u^4-(210/187)*u^2+(175/247)+
					  (70/81)*u^9+(2905/729)*u^7-(35/625482)*u^16-(1085/6561)*u^10+(245/101015343)*
					  u^19+(1295/312741)*u^13-(245/99)*u^8)+(u > 1 & u < 2)*((2870/351)*u^6+(21560/2187)*
					  u^4+(7840/1683)*u^2+(22400/20007)-(15680/6561)*u-(70/81)*u^9+(245/99)*u^8-
					  (11305/2187)*u^7-(980/99)*u^5-(2800/351)*u^3-(245/303046029)*u^19-
					  (665/312741)*u^13+(35/625482)*u^16+(1085/6561)*u^10)}
	   else if (r==1) {fx <- (u <= -1 & u > -2)*((86240/729)*u^2+(28700/117)*u^4+(15680/1683)+
	                     (560/9)*u^7+(10850/729)*u^8+(1400/104247)*u^14+(13720/99)*u^6+(490/1772199)*
						 u^17+(2660/8019)*u^11+(158270/729)*u^5+(19600/99)*u^3+(5600/117)*u)+ 
						 (u <= 0 & u > -1)* ((3920/243)*u^2-(3500/39)*u^4-(420/187)-(560/9)*u^7-
						 (10850/729)*u^8-(1400/104247)*u^14-(13720/99)*u^6-(490/590733)*u^17-
						 (5180/8019)*u^11-(40670/243)*u^5)+(u > 0 & u <= 1)*((3920/243)*u^2-
						 (3500/39)*u^4-(420/187)+(560/9)*u^7-(10850/729)*u^8-(1400/104247)*u^14-
						 (13720/99)*u^6+(490/590733)*u^17+(5180/8019)*u^11+(40670/243)*u^5)+
						 (u > 1 & u < 2)*((86240/729)*u^2+(28700/117)*u^4+(15680/1683)-(560/9)*u^7+
						 (10850/729)*u^8-(158270/729)*u^5+(13720/99)*u^6-(19600/99)*u^3-(5600/117)*u+
						 (1400/104247)*u^14-(2660/8019)*u^11-(490/1772199)*u^17)}
	   else if (r==2) {fx <- (u <= -1 & u > -2)*((114800/39)*u^2+(172480/729)+(3165400/729)*u^3+
	                     (607600/729)*u^6+(26600/729)*u^9+(7840/104247)*u^15+(7840/3)*u^5+
						 (137200/33)*u^4+(39200/33)*u+(19600/8019)*u^12)+ (u <= 0 & u > -1)* 
						 (-(14000/13)*u^2+(7840/243)-(813400/243)*u^3-(607600/729)*u^6-(51800/729)*
						 u^9-(7840/34749)*u^15-(7840/3)*u^5-(137200/33)*u^4-(19600/8019)*u^12)+
						 (u > 0 & u <= 1)*(-(14000/13)*u^2+(7840/243)+(813400/243)*u^3+(51800/729)*
						 u^9-(607600/729)*u^6+(7840/34749)*u^15+(7840/3)*u^5-(137200/33)*u^4-
						 (19600/8019)*u^12)+ (u > 1 & u < 2)*((114800/39)*u^2+(172480/729)-
                         (3165400/729)*u^3+(19600/8019)*u^12-(26600/729)*u^9+(607600/729)*u^6-
						 (7840/3)*u^5-(39200/33)*u+(137200/33)*u^4-(7840/104247)*u^15)}
	   else if (r==3) {fx <- (u <= -1 & u > -2)*( (229600/39)+(212800/81)*u^7+(78400/243)*u^10+
	                     (548800/34749)*u^13+(6076000/243)*u^4+(6330800/243)*u+(156800/3)*u^3+
						 (548800/11)*u^2)+ (u <= 0 & u > -1)* ( -(28000/13)-(414400/81)*u^7-
						 (548800/11583)*u^13-(6076000/243)*u^4-(1626800/81)*u-(78400/243)*u^10-
                         (156800/3)*u^3-(548800/11)*u^2) +(u > 0 & u <= 1)*( (-28000/13)+
						 (414400/81)*u^7+(548800/11583)*u^13+(1626800/81)*u-(6076000/243)*u^4-
						 (78400/243)*u^10+(156800/3)*u^3-(548800/11)*u^2)+ (u > 1 & u < 2)* 
                         ((229600/39)+(78400/243)*u^10-(6330800/243)*u+(6076000/243)*u^4-
						 (548800/34749)*u^13-(212800/81)*u^7-(156800/3)*u^3+(548800/11)*u^2)}
	   else if (r==4) {fx <- (u <= -1 & u > -2)*( -(589568000/81)*u^2-(29478400/33)-(16777600/27)*
	                     u^5-(11603200/3)*u+(2195200/891)*u^11+(784000/27)*u^8-(21952000/3)*u^3-
						 (10976000/3)*u^4)+(u <= 0 & u > -1)* ( (551936000/81)*u^2+(4076800/11)+
						 (2038400/3)*u^5+2822400*u-(2195200/297)*u^11-(784000/27)*u^8+(21952000/3)*
						 u^3+(10976000/3)*u^4)+(u > 0 & u <= 1)*((551936000/81)*u^2+(4076800/11)-
						 (2038400/3)*u^5-2822400*u-(784000/27)*u^8+(2195200/297)*u^11-(21952000/3)*
						 u^3+(10976000/3)*u^4)+ (u > 1 & u < 2)*(-(589568000/81)*u^2-(29478400/33)+
                         (16777600/27)*u^5+(11603200/3)*u+(784000/27)*u^8-(2195200/891)*u^11-
						 (10976000/3)*u^4+(21952000/3)*u^3)}
	   else if (r==5) {fx <- (u <= -1 & u > -2)*((2885120000/81)+(4406080000/27)*u^3+219520000*u^2+
	                     (43904000/27)*u^6+43904000*u^4+137984000*u+(21952000/81)*u^9)+ 
						 (u <= 0 & u > -1)* ( -(1944320000/81)-(486080000/3)*u^3-219520000*u^2-
						 (43904000/27)*u^6-43904000*u^4-125440000*u-(21952000/27)*u^9)+ 
						 (u > 0 & u <= 1)*( -(1944320000/81)+(486080000/3)*u^3-219520000*u^2-
						 (43904000/27)*u^6-43904000*u^4+125440000*u+(21952000/27)*u^9)+
						 (u > 1 & u < 2)*( (2885120000/81)+43904000*u^4+219520000*u^2-137984000*u-
						 (4406080000/27)*u^3-(21952000/81)*u^9+(43904000/27)*u^6)}
	   else if (r==6) {fx <- (u <= -1 & u > -2)*( (-2320640000/3)-(3512320000/3)*u^3-(22798720000/9)*u+
	                     (175616000/9)*u^7-2985472000*u^2+(439040000/9)*u^4)+ (u <= 0 & u > -1)* 
						 ( 689920000+(3512320000/3)*u^3+2540160000*u-(175616000/3)*u^7+2985472000*
						 u^2-(439040000/9)*u^4)+ (u > 0 & u <= 1)*( 689920000-(3512320000/3)*u^3-
						 2540160000*u+(175616000/3)*u^7+2985472000*u^2-(439040000/9)*u^4)+
						 (u > 1 & u < 2)*((-2320640000/3)+(3512320000/3)*u^3+(22798720000/9)*u-
						 (175616000/9)*u^7-2985472000*u^2+(439040000/9)*u^4)}
	   else if (r==7) {fx <- (u <= -1 & u > -2)*( 9834496000+24586240000*u+(2458624000/3)*u^5+
	                     (49172480000/3)*u^2)+(u <= 0 & u > -1)* ( -9834496000-24586240000*u-
						 2458624000*u^5-(49172480000/3)*u^2)+(u > 0 & u <= 1)*( -9834496000+
						 24586240000*u+2458624000*u^5-(49172480000/3)*u^2)+ (u > 1 & u < 2)*
						 (9834496000-24586240000*u-(2458624000/3)*u^5+(49172480000/3)*u^2)}
	   else if (r==8) {fx <- (u <= -1 & u > -2)*(-98344960000*u+(49172480000/3)*u^3-(196689920000/3))+
                  	     (u <= 0 & u > -1)*(98344960000*u-49172480000*u^3+(196689920000/3))+
						 (u > 0 & u <= 1)* (49172480000*u^3-98344960000*u+(196689920000/3))+ 
						 (u > 1 & u < 2)*( 98344960000*u-(49172480000/3)*u^3-(196689920000/3));return(fx)}
	   else if (r==9) {fx <- (u <= -1 & u > -2)*( 98344960000*u+196689920000)+ (u <= 0 & u > -1)*
                    	 (-295034880000*u-196689920000)+(u > 0 & u <= 1)* ( 295034880000*u-
						 196689920000)+ (u > 1 & u < 2)*( -98344960000*u+196689920000)}	  
   }
   else if (kernel=="silverman"){
   r <- r%%4
       if (r==0) {fx <- (u < 0)*((1/16)*exp((1/2)*u*sqrt(2))*(3*sqrt(2)*cos((1/2)*u*sqrt(2))-
	                 3*sqrt(2)*sin((1/2)*u*sqrt(2))+2*sin((1/2)*u*sqrt(2))*u))+ (u >= 0)*((1/16)*
				     exp(-(1/2)*u*sqrt(2))*(3*sqrt(2)*cos((1/2)*u*sqrt(2))+3*sqrt(2)*sin((1/2)*u*
				     sqrt(2))+2*sin((1/2)*u*sqrt(2))*u))}
	   else if (r==1) {fx <- (u < 0)*(-(1/16)*exp((1/2)*u*sqrt(2))*(sqrt(2)*cos((1/2)*u*sqrt(2))+
	                     sqrt(2)*sin((1/2)*u*sqrt(2))-2*cos((1/2)*u*sqrt(2))*u))+ (u >= 0)*
						 (-(1/16)*exp(-(1/2)*u*sqrt(2))*(sqrt(2)*cos((1/2)*u*sqrt(2))+2*cos((1/2)*
						 u*sqrt(2))*u-sqrt(2)*sin((1/2)*u*sqrt(2))))}
	   else if (r==2) {fx <- (u < 0)*((1/16)*exp((1/2)*u*sqrt(2))*(sqrt(2)*cos((1/2)*u*sqrt(2))-
	                     sqrt(2)*sin((1/2)*u*sqrt(2))-2*sin((1/2)*u*sqrt(2))*u))+ (u >= 0)*
						 ((1/16)*exp(-(1/2)*u*sqrt(2))*(sqrt(2)*cos((1/2)*u*sqrt(2))+sqrt(2)*
						 sin((1/2)*u*sqrt(2))-2*sin((1/2)*u*sqrt(2))*u))}
	   else if (r==3) {fx <- (u < 0)*(-(1/16)*exp((1/2)*u*sqrt(2))*(3*sqrt(2)*cos((1/2)*u*sqrt(2))+
	                     3*sqrt(2)*sin((1/2)*u*sqrt(2))+2*cos((1/2)*u*sqrt(2))*u))+ (u >= 0)*
						 (-(1/16)*exp(-(1/2)*u*sqrt(2))*(3*sqrt(2)*cos((1/2)*u*sqrt(2))-3*sqrt(2)*
						 sin((1/2)*u*sqrt(2))-2*cos((1/2)*u*sqrt(2))*u))}
   }
   else if (kernel=="gaussian"){
       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)   
}

####
#### Kernel distribution H(x)
	
			
# kernel_fun_dist <- function(kernel,u)
                         # {
# if (kernel=="epanechnikov")      {H <- function(X) (3/4)*X-(1/4)*X^3+0.5 }
# else if (kernel=="gaussian")     {H <- pnorm}
# else if (kernel=="silverman")    {H <- function(X) (0.25*sqrt(2)*exp(0.5*sqrt(2)*X)*(cos(0.25*pi-0.5*sqrt(2)*X)+sin(0.25*pi-0.5*sqrt(2)*X))) *(X < 0)+ 
#                                                    (-0.25*sqrt(2)*exp(-0.5*sqrt(2)*X)*cos(0.25*pi+0.5*sqrt(2)*X)-0.25*sqrt(2)*exp(-0.5*sqrt(2)*X)*sin(0.25*pi+0.5*sqrt(2)*X)+1)*(X >= 0)}
# else if (kernel=="uniform")      {H <- function(X) 0.5*X+0.5 }
# else if (kernel=="triangular")   {H <- function(X) (0.5+X+0.5*X^2) *(X <0) + (0.5+X+0.5*X^2-X^2) *(X >=0)}
# else if (kernel=="tricube")      {H <- function(X) (0.5+(70/81)*X+(7/81)*X^10 + (10/27)*X^7 + (35/54)*X^4 )*(X < 0)+(0.5+(70/81)*X-(7/81)*X^10 + (10/27)*X^7 - (35/54)*X^4 )*(X >= 0)}
# else if (kernel=="triweight")    {H <- function(X) ((35/32) * X -(35/32)* X^3 +(21/32)* X^5 -(5/32)* X^7 + 0.5)}
# else if (kernel=="biweight")     {H <- function(X) ((15/16) * X -(5/8)* X^3+ (3/16)* X^5 + 0.5) }
# else if (kernel=="cosine")       {H <- function(X) (1/2)*sin((1/2)*pi)+(1/2)*sin((1/2)*pi*X) }

# if (kernel=="gaussian"){Fx <- H(u)}else{
# Fx <- u
# I0 <- (u <= -1)
# I1 <- (u >= 1)
# I2 <- (u > -1 & u < 1)
# Fx[I0] <- 0;Fx[I1] <- 1
# U <- Fx[I2]
# Fx[I2] <- H(U)
# }
# return(Fx)
# }

####
#### 1/(n-1) * sum[H((x-x[j])/h)] ; j != i
#### Author: Alejandro Quintela del Rio and Graciela Estevez Perez
#### Alejandro Quintela del Rio <aquintela@udc.es>

# kernel_dist_wit_i <-function(kernel,y,x,bw)
                    # {
    # n <- length(x)
    # AUX <- matrix(0, n, n)
    # result <- matrix(0,n,length(y))
	# for(j in 1:length(y)){ 
    # AUX <- matrix(rep.int(outer(y[j],x,"-"),n),nrow=n,byrow=TRUE)
	# aux <- kernel_fun_dist(kernel, AUX/bw)
	# diag(aux) <- 0
	# result[,j] <- (1/(n-1))*apply(aux,1,sum)
	                     # }
  	# return(result)
# }


####
#### integral function
#### Author: Alejandro Quintela del Rio and Graciela Estevez Perez
#### Alejandro Quintela del Rio <aquintela@udc.es>

# simp_int <- function (x, fx, n.pts = 256, ret = FALSE) 
# {
    # if (class(fx) == "function") 
        # fx = fx(x)
    # n.x = length(x)
    # if (n.x != length(fx)) 
        # stop("Unequal input vector lengths")
    # if (n.pts < 64) 
        # n.pts = 64
    # ap = approx(x, fx, n = 2 * n.pts + 1)
    # h = diff(ap$x)[1]
    # integral = h * (ap$y[2 * (1:n.pts) - 1] + 4 * ap$y[2 * (1:n.pts)] + ap$y[2 * (1:n.pts) + 1])/3
    # invisible(list(value = sum(integral), cdf = list(x = ap$x[2 * 
        # (1:n.pts)], y = cumsum(integral))))
# }


####
#### \int x K(x) H(x) dx
				
A1_kM <- function(kernel)
        {
if (kernel=="gaussian")          {xKr <- 1/(2*sqrt(pi))}
else if (kernel=="silverman")    {xKr <- (5/32)*sqrt(2)}
else if (kernel=="epanechnikov") {xKr <- 9/70 }
else if (kernel=="uniform")      {xKr <- 1/6 }
else if (kernel=="triangular")   {xKr <- 7/60 }
else if (kernel=="triweight")    {xKr <- 245/2574}
else if (kernel=="tricube")      {xKr <- 4291/39366}
else if (kernel=="biweight")     {xKr <- 25/231}
else if (kernel=="cosine")       {xKr <- 1/8}
return(xKr)
}

####
#### \int x^2 K(x) dx

A2_kM <-function(kernel)
{
if (kernel=="gaussian")          {xKr <- 1}
else if (kernel=="silverman")    {xKr <- 0}
else if (kernel=="epanechnikov") {xKr <- 1/5 }
else if (kernel=="uniform")      {xKr <- 1/3 }
else if (kernel=="triangular")   {xKr <- 1/6 }
else if (kernel=="triweight")    {xKr <- 1/9}
else if (kernel=="tricube")      {xKr <- 35/243}
else if (kernel=="biweight")     {xKr <- 1/7}
else if (kernel=="cosine")       {xKr <- (-8+pi^2)/pi^2}
return(xKr)
}

####
#### \int K(x)^2 dx

A3_kM <-function(kernel)
{
if (kernel=="gaussian")          {xKr <- 1/(2*sqrt(pi))}
else if (kernel=="silverman")    {xKr <- (3/16)*sqrt(2)}
else if (kernel=="epanechnikov") {xKr <- 3/5}
else if (kernel=="uniform")      {xKr <- 1/2}
else if (kernel=="triangular")   {xKr <- 2/3}
else if (kernel=="triweight")    {xKr <- 350/429}
else if (kernel=="tricube")      {xKr <- 175/247}
else if (kernel=="biweight")     {xKr <- 5/7}
else if (kernel=="cosine")       {xKr <- (1/16)*pi^2}
return(xKr)
}

####
#### \int k(x,r)^2 dx

A3_kMr <- function(kernel,r)
        {
 if (kernel=="gaussian"){
     xKr <- integrate(function(x) kernel_fun_der(kernel,x,deriv.order=r)^2, -Inf, Inf)$value}
 else if (kernel =="silverman"){
     xKr <- integrate(function(x) kernel_fun_der(kernel,x,deriv.order=r)^2, -100, +100)$value}else{
     xKr <- integrate(function(x) kernel_fun_der(kernel,x,deriv.order=r)^2, -1, +1)$value}
return(xKr)
}

####
#### \int K''(x)^2 dx	

A4_kM <-function(kernel)
{
if (kernel=="gaussian")          {xKr <- 3/(8*sqrt(pi))}
else if (kernel=="silverman")    {xKr <- (1/16)*sqrt(2)}
else if (kernel=="epanechnikov") {xKr <- 9/2}
else if (kernel=="uniform")      {xKr <- 0}
else if (kernel=="triangular")   {xKr <- 0}
else if (kernel=="triweight")    {xKr <- 35}
else if (kernel=="tricube")      {xKr <- 7840/243}
else if (kernel=="biweight")     {xKr <- 22.5}
else if (kernel=="cosine")       {xKr <- (1/256)*pi^6}
return(xKr)
}			


####
#### \int x^4 * K(x) dx	

A5_kM <-function(kernel)
{
if (kernel=="gaussian")          {xKr <- 3}
else if (kernel=="silverman")    {xKr <- -24}
else if (kernel=="epanechnikov") {xKr <- 3/35}
else if (kernel=="uniform")      {xKr <- 1/5}
else if (kernel=="triangular")   {xKr <- 1/15}
else if (kernel=="triweight")    {xKr <- 1/33}
else if (kernel=="tricube")      {xKr <- 1/22}
else if (kernel=="biweight")     {xKr <- 1/21}
else if (kernel=="cosine")       {xKr <- (pi^4-48*pi^2+384)/pi^4}
return(xKr)
}	

####
#### Derivatives

DD <- function(expr,name, order = 1) {
   if(order < 1)  stop("'order' must be >= 1")
   if(order == 1) D(expr,name)
   else DD(D(expr, name), name, order - 1)
                                     }

####
#### 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)
}

Try the kedd package in your browser

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

kedd documentation built on May 2, 2019, 7:32 a.m.