R/kefJG.R

In SNSequate: Standard and Nonstandard Statistical Models and Methods for Test Equating

###########################################################
# This function pertains to the Ake package 
# (W. E. Wansouwe, S. M. Some and C. C. Kokonendji). 
# Jorge Gonzalez advice the manteiner about a bug that has 
# not yet being solved. Thus, we are literally copying here
# the functions from the Ake package, but solving the bug. 
# Thus, SNSequate do not depend on Ake anymore, but it will
# once the bug is solved. 
# All credit for the functions below is given to the authors 
# of the Ake package. 
############################################################

kef <-
  function(x,t,h,type_data=c("discrete","continuous"),ker=c("bino","triang",
                                                            "dirDU","BE","GA","LN", 
                                                            "RIG"),a0=0,a1=1,a=1,c=2){
    ###########################################################################################################
    # INPUTS:
    #   "x" 	: the target.
    #   "t" 	: the single or the grid value where the function is computed.
    #   "h" 	: the bandwidth parameter.
    #   "ker" 	: the kernel:  "bino" binomial,"triang" discrete triangular,"dirDU" Dirac discrete uniform.
    #			       "BE"  extended beta,"GA" gamma,"LN" lognormal,"RIG" reciprocal inverse Gaussian.     	 
    #   "a0" 	: the left bound of the support of the distribution for extended beta kernel. Default value is 0.
    #   "a1" 	: the right bound of the support of the distribution for extended beta kernel. Default value is 1.
    #   "a" 	: The arm is used only for the discrete triangular distribution.
    #   "c" 	: The number of categories in DiracDU kernel and is used only for DiracDU
    # OUTPUT:
    # Returns the discrete associated kernel value at t.
    ###########################################################################################################
    if (missing(type_data))  stop("argument 'type_data' is omitted")
    if ((type_data=="discrete") & (ker=="GA"||ker=="LN"||ker=="BE" ||ker=="RIG")) 
      stop(" Not appropriate kernel for type_data")
    if ((type_data=="continuous") & (ker=="bino"||ker=="triang"||ker=="dirDU")) 
      stop(" Not appropriate kernel for 'type_data'")
    if ((type_data=="discrete") & missing(ker)) ker<-"bino"
    if ((type_data=="continuous") & missing(ker)) ker<-"GA"
    
    dtrg<-function(x,t,h,a){
      
      if (a==0)
      {
        result <- t
        Logic1 <- (t==x) 
        Logic0 <- (t!=x) 
        result[Logic1]=1
        result[Logic0]=0					
        
      }
      
      else
        
      {	
        
        u=0:a;
        u=sum(u^h)			 
        D=(2*a+1)*(a+1)^h -2*u                 
        result <- t
        Logic0 <- ((t>=(x-a)) & (t<=(x+a)))  # support Sx={x-a,...,x+a} support of the distribution
        Logic1 <- ((t<(x-a))|(t>(x+a)))  
        tval <- result[Logic0]
        result[Logic1]=0
        result[Logic0]<-  ((a+1)^h - (abs(tval-x))^h)/D # Discrete Triangular 				
        
      }
      return(result)
    }
    
    
    
    diracDU<-
      function(x,t,h,c)
      {	
        result<-t
        Logic1 <- (t==x) 
        Logic0 <- (t!=x)
        result[Logic1]<-(1-h)
        result[Logic0]<- (h/(c-1)) 
        
        return(result)        	
        
      }
    
    
    
    if(ker=="bino"){	
      result <- t
      Logic0 <- (t <= x+1) # support Sx={0,1,...,x+1}
      Logic1 <- (x+1 < t)
      tval <- result[Logic0]
      result[Logic1]=0
      result[Logic0]<- dbinom(tval,x+1,(x+h)/(x+1))  # The Binomial kernel 
      
      
    }
    
    
    else	if(ker=="triang"){
      result <- dtrg(x,t,h,a)  # The discrete Triangular kernel
      
      
    }
    
    
    else   if(ker=="dirDU")
    {	
      result <- diracDU(x,t,h,c)  # The Dirac Discrete Uniform kernel
      
    }
    
    
    
    if(ker=="BE"){	
      
      result <- t
      Logic0 <- ((a0<=t)&(t<= a1)) # support 
      Logic1 <- ((t<a0)|(a1<t))
      tval <- result[Logic0]
      result[Logic1]=0
      
      result[Logic0]<- ((1/((a1-a0)^(1+h^(-1))*beta(((x-a0)/((a1-a0)*h))+1,((a1-x)/((a1-a0)*h))+1))))*((tval-a0)^((x-a0)/((a1-a0)*h)))*((a1-tval)^((a1-x)/((a1-a0)*h))) 
      
      
      
    }
    
    
    else if(ker=="GA"){
      result <- t
      Logic0 <- (0<=t) # support 
      Logic1 <- (t<0)
      tval <- result[Logic0]
      result[Logic1]=0		
      result[Logic0]<- dgamma(tval,(x/h)+1,1/h)
      
      
    }
    
    else if(ker=="LN"){ 
      result <- t
      Logic0 <- (0<=t) # support 
      Logic1 <- (t<0)
      tval <- result[Logic0]
      result[Logic1]=0
      # result[Logic0]<-  (1/(tval*h*sqrt(2*pi)))*exp((-1/2)*((1/h)*log(tval/x)-h)^2) 				
      result[Logic0]<- dlnorm(tval,meanlog=log(x)+h^2,sdlog=h)
      
      
    }
    else if(ker=="RIG"){ 
      result <- t
      Logic0 <- (0<t) # support 
      Logic1 <- (t<=0)
      tval <- result[Logic0]
      result[Logic1]<- 0
      eps<-sqrt(x^2+x*h) # see Libengue (2013)
      result[Logic0]<- (1/sqrt(2*pi*h*tval))*exp((-eps/(2*h))*((tval/eps) -2+(eps/tval)))	
      
      
    }
    
    return(result)
  }