# R/penalty.R In LINselect: Selection of Linear Estimators

#### Defines functions calcEDkhiEDkhi2EDkhi1

# ---------------------------------------------------------------
#  LINselect R package
#  Copyright INRA 2017
#  INRA, UR1404, Research Unit MaIAGE
#  F78352 Jouy-en-Josas, France.

###########################################################################
#
# calcul de la penalite
#
###########################################################################

EDkhi1 <- function(x,D,N,q) {
# ---------------------------------------------------------------
# Case : L not too  small
# called by calcEDkhi
# INPUT
# x scalar or vector
# D, N : positive integer
# q : scalar
# OUTPUT
# values in x : scalar or vector
# CALL:
# This function is called by 'calcEDkhi'
# ---------------------------------------------------------------
fct1 <- pf(q=x/(D+2), df1=D+2, df2=N, lower.tail = FALSE)
fct2 <- (x/D)*pf(q=(N+2)*x/(N*D), df1=D, df2=N+2, lower.tail = FALSE)
return(q+fct2-fct1)
} #  fin EDkhi1

EDkhi2  <- function(x,D,N,logq) {
# ---------------------------------------------------------------
# Case : L small
# called by calcEDkhi
# INPUT
# x scalar or vector
# D, N : positive integer
# logq : scalar
# OUTPUT
# values in x : scalar or vector
# CALL:
# This function is called by 'calcEDkhi'
# ---------------------------------------------------------------
t1 <- log(2*(2*x+N*D)/(N*(N+2)*x))-lbeta(1+D/2,N/2)
t3 <- (N/2)*log(N/(N+x))+(D/2)*log(x/(N+x))
return(logq-t1-t3)
} #  fin EDkhi2

penalty <- function
### Calculate the penalty function for estimators selection.
##reference<<  See Baraud  et al. 2010
## \url{http://hal.archives-ouvertes.fr/hal-00502156/fr/} \cr

(Delta, ##<<vector with \code{Dmax}+1 components : weights in the penalty function.
n, ##<<integer : number of observatons.
p, ##<<integer : number of variables.
K ##<<scalar : constant in the penalty function.
) {
# SH 19/02/2013 : penalite et penalty different legerement dans les
# entrees. Il faudra corriger ce doublon
# called by VSelect
# VOIR: AB:08/03/2011 a rajouté l'argument p
# SH 19/02/2013 ok
# ---------------------------------------------------------------
# INPUT
# Delta vector of length Dmax+1
# n,p: integer
# K : scalar
# OUTPUT
# vector of length dmax+1

# CALL:
# This function calls 'calcEDkhi'
# ---------------------------------------------------------------
Dmax <- length(Delta)-1
Dm <- 0:Dmax
pen <- Dm*0
Nm <- n-Dm
Lm <- Delta
EDkhi <- calcEDkhi(p, Dm+1,Nm-1,Lm)
if ( EDkhi$err != 0) { # error or warning cases mess <- paste(EDkhi$mess,
"\n err=", EDkhi$err, "n n=", n, "Dm=", Dm) if (EDkhi$err ==2) {
##note<< The values of the penalty function greater than
##1e+08 are set to 1e+08.
warning(mess)
}
else
stop(mess)
##note<<  If for some \code{Delta(d)} the
##equation \eqn{\phi}\code{(x) = exp(-Delta(d)/(d+1))}  has no
##solution, then the execution is stopped.
} # fin err
EDkhiR <- EDkhi$EDkhi2 pen <- K*(Nm/(Nm-1))*EDkhiR names(pen) <- Dm return(pen) ### A vector with the same length as Delta: for each \code{d}=0, ..., \code{Dmax}, let ### \code{N}=\code{n}-\code{d}, \code{D}=\code{d+1} and \cr ### \code{pen(d) = x K N/(N-1)} where \code{x} satisfies ###\cr ###\cr ### \eqn{\phi}\code{(x) = exp(-Delta(d))}, when \code{Delta(d)<50}, ### \cr ### where \eqn{\phi}\code{(x)=pf(q=x/(D+2),df1=D+2,df2=N-1,lower.tail=F)-(x/D)pf(q=(N+1)x/D(N-1),df1=D,df2=N+1,lower.tail=F)} ###\cr ###\cr ### \eqn{\psi}\code{(x) = Delta(d)}, when ### \code{Delta(d)}\eqn{\ge}\code{50}, \cr ### where \eqn{\psi}\code{(x)=lbeta(1+D/2,(N-1)/2)-log(2(2x+(N-1)D)/((N-1)(N+2)x))-(N-1)/2log((N-1)/(N-1+x))-(D/2)log(x/(N-1+x)) ### } } # fin penalty calcEDkhi <- function(p, D,N,L) # VOIR: AB:08/03/2011 a rajouté l'argument p # EDkhi <- calcEDkhi(Dm+1,Nm-1,LmC) # called by penalite and penalty # D<-Dm+1; N=Nm-1; L=Lm # SH : Nouvelle version # --------------------------------------------------------------- # INPUT # p : integer # D, N vectors of integers # K : vector # OUTPUT # array : list with components # EDkhi2 value of EDkhi # err : integer. Error or Warning code. # mess : error message if err not equal 0 # CALL: # This function calls 'EDkhi1', 'EDkhi2'. # It is called by 'penalty' # --------------------------------------------------------------- { Bsup <- 1e+08 err <- 0 mess <- NULL xq <- 0*D for (i in 2:length(D)) { if (xq[i-1] > Bsup) { xq[i:length(D)] <- Inf err <- 2 mess <- paste("the values of the penalty function greater than", Bsup, "are set to Inf") break } if (i==2) { xInf <- 0 if (N[2]>5) { Delta <- (L[2]+log(5)+1/N[2])/(1-5/N[2]) U <- sqrt((1+2*D[2]/(N[2]+2))*2*Delta/D[2]) xSup <- D[2]*(1+exp(2*Delta/(N[2]+2))*U)**2 } if (N[2]<=5) { xSup <- seq(10,100,by=10)*p fxSup <- EDkhi1(xSup,D[i],N[i],exp(-L[i])) if (max(fxSup) <0 ) { err <- 4 mess <- paste("the values of the penalty function cannot be calculated\n", xq) break } if (max(fxSup) >=0 ) { xSup <- min(xSup[fxSup>=0]) } } } else { # (i>2) xInf <- xq[i-1] xSup <- xq[i-1]*seq(10,100,by=10) fxSup <- EDkhi1(xSup,D[i],N[i],exp(-L[i])) if (max(fxSup) <0 ) { err <- 3 mess <- paste("the values of the penalty function cannot be calculated\n", xq) break } xSup <- min(xSup[fxSup>=0]) } if (L[i] < 50) { xx <- try(uniroot(EDkhi1,lower=xInf,upper=xSup,D=D[i], N=N[i],q=exp(-L[i]))) if (!is.list(xx)) { err <- 1 mess <- xx break } xq[i] <- xx$root
}
else {
# L[i] >= 50
fxSup <- EDkhi2(xSup,D[i],N[i],-L[i])
if (fxSup <0) xSup<- 2*xSup
z <-
try(uniroot(EDkhi2,lower=xInf,upper=xSup,D=D[i],
N=N[i],logq=-L[i]))
if (!is.list(xx)) {
err <- 1
mess <- xx
break
}
xq[i] <- z\$root
}
}
return(list(EDkhi2=xq,err=err, mess=mess))
} #  fin calcEDkhi


## Try the LINselect package in your browser

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

LINselect documentation built on May 2, 2019, 4:09 a.m.