Nothing
# ---------------------------------------------------------------
# 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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.