Mspline/ 40755 20207 20154 0 7615722557 11523 5 ustar cantoni metri Mspline/Ins3.4 100744 20207 20154 310 7614241012 12451 0 ustar cantoni metri #! /bin/sh
S=${S-Splus}
mkdir .Data .Data/.Help
cat Msplinesfun3.4.s | $S
for f in help3.4/*.d
do
mv $f .Data/.Help/`basename $f .d`
done
$S help.findsum .Data
rm -f .Data/.Audit .Data/.Last.value
Mspline/help6.0/ 40755 20207 20154 0 7615723304 12666 5 ustar cantoni metri Mspline/help6.0/frob.sgml 100644 20207 20154 10356 7515013360 14615 0 ustar cantoni metri
]
>
frob
Fit an M-type smoothing spline to the input data.
Fit an M-type smoothing spline as described in Eva Cantoni and Elvezio Ronchetti,
Resistant Selection of the Smoothing Paramater for Smooting Splines, 2001,
Statistics and Computing, 11, 141-146.
frob(xx, yy, lambda, chuber=1.345, n.orig=length(xx), compteurmax=500, mytol=.Machine$double.eps^0.25, Smat=F)
Values of the predictor variable.
Response variable, of the same length as xx.
Smoothing parameter corresponding to the formulation of the reference paper.
Note that it is not the same parameter as spar in the function smooth.spline()
in S-PLUS. In fact spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
Tuning constant of the Huber psi function. The default value is set to 1.345.
Number of original observations. This parameter controls the weighting procedure
when there are ties and is used to assess the condition sum(w)=n.orig in
my.smooth.spline(). The default is set to length(xx).
Maximum number of iterations of the fitting algorithm.
Convergence threshold.
Logical value indicating wheter the whole smoother matrix has to be computed.
The default is F. In this case, only the diagonal elements of the smoother matrix
are computed. Moreover, the RCp criterion is not computed, and only the
cross-validation score is given.
A list is returned with the following components:
The fitted M-type smoothing spline corresponding to the ordered xx.
The associated estimation of scale (Huber Proposal 2).
The external estimation of scale used in the construction of RCp (if Smat=T).
The ordered distinct values of xx.
The y-values used at the unique x values (weighted averages of the input yy).
Value of match(xx, unique(sort(xx))). This allow one to recover the vector of
fitted values with respect to the original xx by considering estim[myo].
weights used in the fit. This has the same length as xx, and in the case of ties,
will consist of the accumulated weights at each unique value of x.
The smoother matrix S, if Smat=T.
The diagonal elements of the smoother matrix S.
The expectation of the derivative of the psi function.
The expectation of the squared psi function.
The value of the robust cross-validation criterion.
The value of the robust Cp criterion, if Smat=T.
The tuning constant used for the fit.
The smoothing parameter used for the fit.
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
opt.cv, opt.RCp.
myfit<- frob(myx,myy,lambda=0.01)
motif()
plot(myx,myy)
lines(myfit$x,myfit$estim)
robust,
function
Mspline/help6.0/matS.sgml 100644 20207 20154 2014 7515013360 14541 0 ustar cantoni metri
]
>
matS
Auxiliary function for frob().
Auxiliary function for frob() which computes the smoother matrix S.
matS(xx, lambda, wei, n.orig)
Values of the predictor variable.
Smoothing parameter.
weights used in the fit.
Number of original observations.
The smoother matrix S.
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
frob.
function
Mspline/help6.0/opt.RCp.sgml 100644 20207 20154 4275 7515013360 15135 0 ustar cantoni metri
]
>
opt.RCp
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust Cp.
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
Cp as described in Eva Cantoni and Elvezio Ronchetti, Resistant Selection of
the Smoothing Paramater for Smooting Splines, 2001, Statistics and Computing,
11, 141-146.
opt.RCp(xx, yy, bornes.opt, ind.chuber=1.345)
Values of the predictor variable.
Response variable, of the same length as xx.
Interval in which the optimal smoothing parameter has to be found.
Tuning constant of the Huber psi function. The default value is set to 1.345.
A list is returned with the following components:
The optimal value of the smoothing parameter obtained by minimization of the
robust Cp criterion.
The value of the robust Cp for the optimal value of the smoothing parameter.
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
frob,opt.cv.
myoptimalfit<- opt.RCp(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
robust,
function
Mspline/help6.0/opt.cv.sgml 100644 20207 20154 4365 7515013360 15061 0 ustar cantoni metri
]
>
opt.cv
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust cross-validation.
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
cross-validation as described in Eva Cantoni and Elvezio Ronchetti, Resistant
Selection of the Smoothing Paramater for Smooting Splines, 2001, Statistics and
Computing, 11, 141-146.
opt.cv(xx, yy, bornes.opt, ind.chuber=1.345)
Values of the predictor variable.
Response variable, of the same length as xx.
Interval in which the optimal smoothing parameter has to be found.
Tuning constant of the Huber psi function. The default value is set to 1.345.
A list is returned with the following components:
The optimal value of the smoothing parameter obtained by minimization of the
robust cross-validation criterion.
The value of the robust cross-validation for the optimal value of the smoothing
parameter.
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
frob,opt.RCp.
myoptimalfit<- opt.cv(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
robust,
function
Mspline/help6.0/psihuber.sgml 100644 20207 20154 1171 7515013360 15461 0 ustar cantoni metri
]
>
psihuber
Huber psi function.
Define a Huber psi function.
psihuber(x, chuber)
Value at which the function needs to be evaluated.
Value of the tuning constant.
Value of the function at point x.
function
Mspline/help6.0/e.psi.sgml 100644 20207 20154 2465 7615723304 14674 0 ustar cantoni metri
]
>
e.psi
Auxiliary function for frob().
Auxiliary function for frob(), which computes the expectations terms (at the
normal model) appearing in the formulation of an M-type smoothing spline and in
the resistant criterion (RCV and RCp) for the selection of the smoothing parameter
e.psi(chuber=1.345)
the tuning constant of the Huber psi function.
The function returns the expectations (at the normal model). See Eva Cantoni and
Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting
Splines", Statistics and Computing, 11, 141-146, for more details.
The integrations are performed numerically via integrate().
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
frob, opt.cv, opt.RCp, integrate.
function
Mspline/help6.0/my.smooth.spline.sgml 100644 20207 20154 3357 7515013360 17076 0 ustar cantoni metri
]
>
my.smooth.spline
Auxiliary function for frob().
Auxiliary function for frob(). It is a modified version of smooth.spline, which
assumes that ties are treated before. (This is done in the first part of frob()).
my.smooth.spline(x, y, w, n.orig, df=5, spar=0, cv=F, all.knots=F, df.offset=0,
penalty=1)
Values of the predictor variable.
Response variable, of the same length as x.
Weights used in the fit.
Number of original observations.
Smoothing parameter satisfying spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
The arguments of my.smooth.spline listed below are in general not used. They are
present because my.smooth.spline is a modified version of smooth.spline.
This function returns an object of class smooth.spline. See the help of this
function for further details.
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
frob.
function
Mspline/helpHTML/ 40755 20207 20154 0 7615730205 13125 5 ustar cantoni metri Mspline/helpHTML/e.psi.html 100755 20207 20154 4070 7615724127 15140 0 ustar cantoni metri
e.psi
Auxiliary function for frob().
DESCRIPTION:
Auxiliary function for frob(), which computes the expectations terms (at the
normal model) appearing in the formulation of an M-type smoothing spline and in
the resistant criterion (RCV and RCp) for the selection of the smoothing parameter
USAGE:
e.psi(chuber=1.345)
REQUIRED ARGUMENTS:
- chuber
-
the tuning constant of the Huber psi function.
VALUE:
The function returns the expectations (at the normal model). See Eva Cantoni and
Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting
Splines", 2001, Statistics and Computing, 11, 141-146, for more details.
DETAILS:
The integrations are performed numerically via integrate().
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
SEE ALSO:
frob, opt.cv, opt.RCp, integrate.
Mspline/helpHTML/frob.html 100755 20207 20154 11374 7615725337 15103 0 ustar cantoni metri
frob
Fit an M-type smoothing spline to the input data.
DESCRIPTION:
USAGE:
frob(xx, yy, lambda, chuber=1.345, n.orig=length(xx), compteurmax=500., mytol=.Machine$double.eps^0.25, Smat=F)
REQUIRED ARGUMENTS:
- xx
-
Values of the predictor variable.
- yy
-
Response variable, of the same length as xx.
- lambda
-
Smoothing parameter corresponding to the formulation of the reference paper.
Note that it is not the same parameter as spar in the function smooth.spline()
in S-PLUS. In fact spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
OPTIONAL ARGUMENTS:
- chuber
-
Tuning constant of the Huber psi function. The default value is set to 1.345.
- n.orig
-
Number of original observations. This parameter controls the weighting procedure
when there are ties and is used to assess the condition sum(w)=n.orig in
my.smooth.spline(). The default is set to length(xx).
- compteurmax
-
Maximum number of iterations of the fitting algorithm.
- mytol
-
Convergence threshold.
- Smat
-
Logical value indicating wheter the whole smoother matrix has to be computed.
The default is F. In this case, only the diagonal elements of the smoother matrix
are computed. Moreover, the RCp criterion is not computed, and only the
cross-validation score is given.
VALUE:
A list is returned with the following components:
estim
The fitted M-type smoothing spline corresponding to the ordered xx.
sigmahat
The associated estimation of scale (Huber Proposal 2).
sigma.ext
The external estimation of scale used in the construction of RCp (if Smat=T).
x
The ordered distinct values of xx.
yin
The y-values used at the unique x values (weighted averages of the input yy).
myo
Value of match(xx, unique(sort(xx))). This allow one to recover the vector of
fitted values with respect to the original xx by considering estim[myo].
weights
weights used in the fit. This has the same length as xx, and in the case of ties,
will consist of the accumulated weights at each unique value of x.
Smatrix
The smoother matrix S, if Smat=T.
diagS
The diagonal elements of the smoother matrix S.
epsiprime
The expectation of the derivative of the psi function.
epsicarre
The expectation of the squared psi function.
cv.score
The value of the robust cross-validation criterion.
RCp
The value of the robust Cp criterion, if Smat=T.
chuber
The tuning constant used for the fit.
lambda
The smoothing parameter used for the fit.
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
SEE ALSO:
opt.cv, opt.RCp.
EXAMPLES:
myfit<- frob(myx,myy,lambda=0.01)
motif()
plot(myx,myy)
lines(myfit$x,myfit$estim)
Mspline/helpHTML/matS.html 100755 20207 20154 3365 7615726011 15026 0 ustar cantoni metri
matS
Auxiliary function for frob().
DESCRIPTION:
Auxiliary function for frob() which computes the smoother matrix S.
USAGE:
matS(xx, lambda, wei, n.orig)
REQUIRED ARGUMENTS:
- xx
-
Values of the predictor variable.
- lambda
-
Smoothing parameter.
- wei
-
weights used in the fit.
- n.orig
-
Number of original observations.
VALUE:
The smoother matrix S.
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
SEE ALSO:
frob.
Mspline/helpHTML/my.smooth.spline.html 100755 20207 20154 5062 7615726507 17356 0 ustar cantoni metri
my.smooth.spline
DESCRIPTION:
USAGE:
my.smooth.spline(x, y, w, n.orig, df=5., spar=0., cv=F, all.knots=F, df.offset=0., penalty=1.)
REQUIRED ARGUMENTS:
- x
-
Values of the predictor variable.
- y
-
Response variable, of the same length as x.
- w
-
Weights used in the fit.
- n.orig
-
Number of original observations.
- spar
-
Smoothing parameter satisfying spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
OPTIONAL ARGUMENTS:
The arguments of my.smooth.spline listed below are in general not used. They are
present because my.smooth.spline is a modified version of smooth.spline.
- df
-
- cv
-
- all.knots
-
- df.offset
-
- penalty
-
VALUE:
This function returns an object of class smooth.spline. See the help of this
function for further details.
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
Mspline/helpHTML/opt.RCp.html 100755 20207 20154 5605 7615727506 15420 0 ustar cantoni metri
opt.RCp
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust Cp.
DESCRIPTION:
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
Cp as described in Eva Cantoni and Elvezio Ronchetti, Resistant Selection of
the Smoothing Paramater for Smooting Splines, 2001, Statistics and Computing,
11, 141-146.
USAGE:
opt.RCp(xx, yy, bornes.opt, ind.chuber=1.345)
REQUIRED ARGUMENTS:
- xx
-
Values of the predictor variable.
- yy
-
Response variable, of the same length as xx.
- bornes.opt
-
Interval in which the optimal smoothing parameter has to be found.
OPTIONAL ARGUMENTS:
- ind.chuber
-
Tuning constant of the Huber psi function. The default value is set to 1.345.
VALUE:
A list is returned with the following components:
lambda.opt
The optimal value of the smoothing parameter obtained by minimization of the
robust Cp criterion.
cv.opt
The value of the robust Cp for the optimal value of the smoothing parameter.
ajust.opt
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
SEE ALSO:
frob,opt.cv
EXAMPLES:
myoptimalfit<- opt.RCp(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
Mspline/helpHTML/opt.cv.html 100755 20207 20154 5702 7615730037 15333 0 ustar cantoni metri
opt.cv
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust cross-validation.
DESCRIPTION:
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
cross-validation as described in Eva Cantoni and Elvezio Ronchetti, Resistant
Selection of the Smoothing Paramater for Smooting Splines, 2001, Statistics and
Computing, 11, 141-146.
USAGE:
opt.cv(xx, yy, bornes.opt=c(0., 1.), ind.chuber=1.345)
REQUIRED ARGUMENTS:
- xx
-
Values of the predictor variable.
- yy
Response variable, of the same length as xx.
-
- bornes.opt
-
Interval in which the optimal smoothing parameter has to be found.
OPTIONAL ARGUMENTS:
- ind.chuber
-
Tuning constant of the Huber psi function. The default value is set to 1.345.
VALUE:
A list is returned with the following components:
lambda.opt
The optimal value of the smoothing parameter obtained by minimization of the
robust cross-validation criterion.
cv.opt
The value of the robust cross-validation for the optimal value of the smoothing parameter.
ajust.opt
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
REFERENCES:
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing
Paramater for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
SEE ALSO:
frob,opt.RCp
EXAMPLES:
myoptimalfit<- opt.cv(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
Mspline/helpHTML/psihuber.html 100755 20207 20154 2235 7615730205 15736 0 ustar cantoni metri
psihuber
Huber psi function.
DESCRIPTION:
Define a Huber psi function.
USAGE:
psihuber(x, chuber)
REQUIRED ARGUMENTS:
- x
-
Value at which the function needs to be evaluated.
- chuber
-
Value of the tuning constant.
VALUE:
Value of the function at point x.
Mspline/Ins6.0 100744 20207 20154 216 7614241011 12454 0 ustar cantoni metri #! /bin/sh
S=${S-Splus}
$S CHAPTER
cat Msplinesfun6.0.s | $S
cd help6.0
for f in *.sgml
do
$S HINSTALL ../.Data $f
done
cd ..
$S BUILD_JHELP
Mspline/Msplinesfun3.4.s 100644 20207 20154 36144 7614237644 14561 0 ustar cantoni metri e.psi_function(chuber=1.345)
{
# This function computes by numerical integration the expectations terms
# appearing in the formulation of an M-type smoothing spline and in the
# resistant criterion (RCV and RCp) for the selection of the smoothing parameter
assign("chuber",chuber,frame=1)
basepsi_function(x)
{
x*pmin(1,chuber/abs(x))
}
assign("basepsi",basepsi,frame=1)
basepsiprime_function(x)
{
1*(abs(x)compteurmax)
stop(paste("Convergence not obtained after",compteurmax,"iterations"))
}
# Cross-validation score. We suppose psi(0)=0.
cv.score_sigmahat.new^2/epsiprime^2*weighted.mean((psihuber((yy-frob.new[myo])
/sigmahat.new,chuber)/(1-(diag.matsol[myo]*rep(1,length(yy)))/ww[myo]))^2,
rep(1,length(yy)))
sigma.ext_NULL
RCp_NULL
matsol_NULL
if(Smat)
{
Ident_diag(1,nrow=length(sx))
matsol_matS(sx,lambda*sigmahat.old/epsiprime,wei=ww,n.orig)
# External estimation of scale for the RCp criterion
acoeff_diff(sx)[-1]/diff(sx,lag=2)
bcoeff_diff(sx)[-(length(sx)-1)]/diff(sx,lag=2)
normcoeff_sqrt(ww[-c(length(sy)-1,length(sy))]^2*acoeff^2+
ww[-c(1,2)]^2*bcoeff^2+ww[-c(1,length(sy))]^2)
sigma.ext_scale.tau((acoeff*ww[-c(length(sy)-1,length(sy))]*
sy[-c(length(sy)-1,length(sy))]+bcoeff*ww[-c(1,2)]*sy[-c(1,2)]-
ww[-c(1,length(sy))]*sy[-c(1,length(sy))])/normcoeff,center=0)
# RCp criterion
RCp.RASR_sum(psihuber((yy-frob.new[myo])/sigma.ext,chuber)^2)
UmoinsV_nb*epsicarre - 2/epsiprime*epsiprimepsicarre*sum(diag(matsol)) +
1/sigma.ext^2*(epsiprimecarre - epsicarresurxcarre)*
sum(diag((matsol-Ident)%*%frob.new%*%t(frob.new)%*%t(matsol-Ident)))+
1/epsiprime^2*(epsicarrepsiprimecarre - epsiquatsurxcarre)*
sum(diag(matsol%*%t(matsol)))
RCp_RCp.RASR-UmoinsV
}
# Output
invisible(return(estim=frob.new,sigmahat=sigmahat.new,sigma.ext=sigma.ext,x=sx,
yin=sy,myo=myo,weights=ww,Smatrix=matsol,diagS=diag.matsol,epsiprime=epsiprime,
epsicarre=epsicarre,cv.score=cv.score,RCp=RCp,chuber=chuber,lambda=lambda))
}
matS_function(xx,lambda,wei,n.orig)
{
# Computation of the matrix defining a resistant spline.
nb_length(xx)
matS_matrix(ncol=nb,nrow=nb)
Ident_diag(1,nrow=nb)
myspar_lambda/(max(xx)-min(xx))^3/mean(wei)
for(j in 1:nb)
{
matS[,j]_my.smooth.spline(xx,Ident[,j],w=wei,n.orig=n.orig,spar=myspar,all.knots=F)$y
}
matS
}
my.smooth.spline_function(x, y, w , n.orig , df = 5, spar = 0, cv = F,
all.knots = F, df.offset = 0, penalty = 1)
#This function assumes that ties are treated before.
{
my.call <- match.call()
if(missing(y)) {
if(is.list(x)) {
if(any(is.na(match(c("x", "y"), names(x)))))
stop("cannot find x and y in list")
y <- x$y
x <- x$x
}
else if(is.complex(x)) {
y <- Im(x)
x <- Re(x)
}
else if(is.matrix(x) && ncol(x) == 2) {
y <- x[, 2]
x <- x[, 1]
}
else {
y <- x
x <- time(x)
}
}
n <- length(x)
if(n != length(y))
stop("lengths of x and y must match")
if(missing(w))
w <- rep(1, n)
else {
if(length(x) != length(w))
stop("lengths of x and w must match")
if(any(w < 0))
stop("all weights should be non-negative")
#w <- (w * sum(w > 0))/sum(w) #this makes the weights sum to n
# Eva: This line was replaced by the following one, in order to construct
# the matrix S correctly, when there are ties in the x's (with matS.s)
if (sum(w)!=n.orig) stop("sum of weigths muste be equal to n.orig")
}
if(missing(spar))
ispar <- 0
else {
if(spar < 1.01e-15)
ispar <- 0
else ispar <- 1
}
if(cv) {
icrit <- 2
}
else {
icrit <- 1
}
dfinfo <- df.offset
if(!missing(df)) {
if(df > 1 & df < n) {
icrit <- 3
dfinfo <- df
}
else warning("you must supply 1 < df < n")
}
n <- as.integer(n)
isetup <- as.integer(0)
ier <- as.integer(1)
mode(x) <- "double"
mode(y) <- "double"
mode(w) <- "double"
mode(spar) <- "double"
mode(ispar) <- "integer"
mode(icrit) <- "integer"
mode(isetup) <- "integer"
mode(ier) <- "integer"
x <- signif(x, 6)
sx <- unique(sort(x))
o <- match(x, sx)
nx <- length(sx)
# Eva: this part is already done in frob.s
# collaps <- .Fortran("suff",
# n,
# as.integer(nx),
# as.integer(o),
# x,
# y,
# w,
# double(nx),
# ybar = double(nx),
# wbar = double(nx),
# double(n))
xbar <- (sx - sx[1])/(sx[nx] - sx[1])
if(all.knots) {
knot <- c(rep(xbar[1], 3), xbar, rep(xbar[nx], 3))
nk <- nx + 2
}
else {
knot <- .Fortran("sknotl",
xbar,
as.integer(nx),
knot = double(n + 6),
k = integer(1))
nk <- knot$k - 4
knot <- knot$knot[seq(knot$k)]
}
low.parm <- 0
high.parm <- 1.5
fit <- .Fortran("qsbart",
as.double(penalty),
as.double(dfinfo),
x = xbar,
y = y, # y = collaps$ybar
w = w, # w = collaps$wbar
as.integer(nx),
knot,
as.integer(nk),
coef = double(nk),
ty = double(nx),
lev = double(nx),
crit = double(1),
c(icrit, ispar),
spar = spar,
c(as.double(low.parm), as.double(high.parm), as.double(0.001)),
isetup,
double((17 + nk) * nk),
as.integer(4),
as.integer(1),
ier = ier)
if(fit$ier > 0) {
warning("smoothing parameter value too small or too large")
fit$ty <- rep(mean(y), nx)
}
lev <- fit$lev
df <- sum(lev)
win <- w[1:nx] # win <- collaps$wbar[1:nx]
if(cv) {
ww <- win
ww[!(ww > 0)] <- 1 # for the next line
cv.crit <- weighted.mean(((y - fit$ty[o])/(1 - (lev[o] * w)/ww[
o]))^2, w)
}
else cv.crit <- weighted.mean((y - fit$ty[o])^2, w)/(1 - (df.offset +
penalty * df)/sum(win))^2
# pen.crit <- sum(win * (collaps$ybar - fit$ty) * collaps$ybar)
pen.crit <- sum(win * (y - fit$ty) * y)
fit.object <- list(knot = knot, nk = nk, min = sx[1], range = sx[nx] -
sx[1], coef = fit$coef)
class(fit.object) <- "smooth.spline.fit"
object <- list(x = sx, y = fit$ty, w = win, yin = y[1:nx],
#yin = collaps$ybar[1:nx]
lev = lev, cv.crit = cv.crit, pen.crit = pen.crit, df = df,
spar = fit$spar, fit = fit.object, call = my.call)
class(object) <- "smooth.spline"
object
}
opt.RCp_function(xx,yy,bornes.opt,ind.chuber=1.345)
{
# Fit of an M-type smoothing spline with parameter chosen by robust C_p,
# as described in Eva Cantoni, Elvezio
# Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting Splines",
# Cahier du Département d'Econométrie, University of Geneva, 98.06, August 1998
if(length(xx) != length(yy))
stop("lengths of x and y must match")
nb_length(xx)
# Expectation terms for the RCp criterion
psi.integr_e.psi(ind.chuber)
epsiprime_psi.integr$epsiprime
epsicarre_psi.integr$epsicarre
epsicarresurxcarre_psi.integr$epsicarresurxcarre
epsiprimepsicarre_psi.integr$epsiprimepsicarre
epsiquatsurxcarre_psi.integr$epsiquatsurxcarre
epsiprimecarre_psi.integr$epsiprimecarre
epsicarrepsiprimecarre_psi.integr$epsicarrepsiprimecarre
assign("ind.chuber",ind.chuber,frame=1)
assign("bornes.opt",bornes.opt,frame=1)
sx_unique(sort(xx))
myo_match(xx,sx)
Ident_diag(1,length(sx))
assign("Ident",Ident,frame=1)
# Objective function to minimize.
RCp.value_function(lambda,xx,yy,nb,epsiprime,epsicarre,epsiprimepsicarre,
epsiprimecarre,epsicarresurxcarre,epsicarrepsiprimecarre,epsiquatsurxcarre,myo)
{
ajust_frob(xx,yy,lambda,chuber=ind.chuber,Smat=T)
fchap_ajust$estim
Smatrix_ajust$Smatrix
sigma.ext_ajust$sigma.ext
ww_ajust$weights
RCp.RASR_sum(psihuber((yy-fchap[myo])/sigma.ext,ind.chuber)^2)
UmoinsV_nb*epsicarre - 2/epsiprime*epsiprimepsicarre*sum(diag(Smatrix)) +
1/sigma.ext^2*(epsiprimecarre - epsicarresurxcarre)*
sum(diag((Smatrix-Ident)%*%fchap%*%t(fchap)%*%t(Smatrix-Ident)))+
1/epsiprime^2*(epsicarrepsiprimecarre - epsiquatsurxcarre)*
sum(diag(Smatrix%*%t(Smatrix)))
RCp.RASR-UmoinsV
}
# Minimisation
opt.lambda_optimize(RCp.value,bornes.opt,xx=xx,yy=yy,nb=nb,epsiprime=epsiprime,
epsicarre=epsicarre,epsiprimepsicarre=epsiprimepsicarre,
epsiprimecarre=epsiprimecarre,
epsicarresurxcarre=epsicarresurxcarre,
epsicarrepsiprimecarre=epsicarrepsiprimecarre,
epsiquatsurxcarre=epsiquatsurxcarre,myo=myo)
# Optimal values.
lambda.opt_opt.lambda$minimum
RCp.opt_opt.lambda$objective
# Checks
mytol <- .Machine$double.eps^.25
if (bornes.opt[2]-lambda.opt < mytol)
warning(paste("lambda.opt=",round(lambda.opt,6),"too close to upper bound",
bornes.opt[2], "\n Search interval need to be widen"))
if(ind.chuber>10000 && lambda.opt < mytol)
stop("lambda.opt cannot be obtained due to the instability of Cp. \n See comment of Figure 3, p. 145, in Cantoni & Ronchetti (2001)")
if (lambda.opt - bornes.opt[1] < mytol)
warning(paste("lambda.opt=",round(lambda.opt,6),"too close to lower bound",
bornes.opt[1], "\n Search interval need to be widen"))
# Fit for the optimal value of the smoothing parameter
ajust.opt_frob(xx,yy,lambda.opt,chuber=ind.chuber,Smat=T)
# Output
invisible(return(lambda.opt,RCp.opt,ajust.opt))
}
opt.cv_function(xx,yy,bornes.opt=c(0,1),ind.chuber=1.345)
{
# Fit of an M-type smoothing spline with parameter chosen by robust
# cross-validation, as described in Eva Cantoni, Elvezio
# Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting Splines",
# Cahier du Département d'Econométrie, University of Geneva, 98.06, August 1998
if(length(xx) != length(yy))
stop("lengths of x and y must match")
nb <- length(xx)
psi.integr_e.psi(ind.chuber)
epsiprime_psi.integr$epsiprime
assign("ind.chuber",ind.chuber,frame=1)
assign("bornes.opt",bornes.opt,frame=1)
# Distinct values of xx.
sx_unique(sort(xx))
myo_match(xx,sx)
# Objective function to minimize.
cv.value_function(lambda,xx,yy,nb,epsiprime,myo)
{
ajust_frob(xx,yy,lambda,chuber=ind.chuber)
fchap_ajust$estim
diagS_ajust$diagS
sigmahat_ajust$sigmahat
ww_ajust$weights
sigmahat^2/epsiprime^2*weighted.mean((psihuber((yy-fchap[myo])/sigmahat,
ind.chuber)/(1-(diagS[myo]*rep(1,length(yy)))/ww[myo]))^2,rep(1,length(yy)))
}
# Minimisation
opt.lambda_optimize(cv.value,bornes.opt,xx=xx,yy=yy,nb=nb,epsiprime=epsiprime,
myo=myo)
# Optimal values.
lambda.opt_opt.lambda$minimum
cv.opt_opt.lambda$objective
# Checks
mytol <- .Machine$double.eps^.25
if (bornes.opt[2]-lambda.opt < mytol)
warning(paste("lambda.opt=",round(lambda.opt,6),"too close to upper bound",
bornes.opt[2], "\n Search interval need to be widen"))
if (lambda.opt - bornes.opt[1] < mytol)
warning(paste("lambda.opt=",round(lambda.opt,6),"too close to lower bound",
bornes.opt[1], "\n Search interval need to be widen"))
# Fit for the optimal value of the smoothing parameter
ajust.opt_frob(xx,yy,lambda.opt,chuber=ind.chuber)
# Output
invisible(return(lambda.opt,cv.opt,ajust.opt))
}
psihuber_function(x,chuber)
{
# This function defines the psi Huber function,
# with tuning constant equal to chuber.
x*pmin(1,chuber/abs(x))
}
Mspline/Msplinesfun6.0.s 100644 20207 20154 37041 7614253460 14547 0 ustar cantoni metri e.psi <- function(chuber = 1.345)
{
# This function computes by numerical integration the expectations terms
# appearing in the formulation of an M-type smoothing spline and in the
# resistant criterion (RCV and RCp) for the selection of the smoothing parameter
assign("chuber", chuber, frame = 1.)
basepsi <- function(x)
{
x * pmin(1., chuber/abs(x))
}
assign("basepsi", basepsi, frame = 1.)
basepsiprime <- function(x)
{
1. * (abs(x) < chuber)
}
assign("basepsiprime", basepsiprime, frame = 1.)
psiprime <- function(x)
{
basepsiprime(x) * dnorm(x)
}
assign("psiprime", psiprime, frame = 1.)
psiprimecarre <- function(x)
{
basepsiprime(x) * basepsiprime(x) * dnorm(x)
}
assign("psiprimecarre", psiprimecarre, frame = 1.)
psicarre <- function(x)
{
basepsi(x) * basepsi(x) * dnorm(x)
}
assign("psicarre", psicarre, frame = 1.)
psiprimepsicarre <- function(x)
{
basepsi(x) * basepsi(x) * basepsiprime(x) * dnorm(x)
}
assign("psiprimepsicarre", psiprimepsicarre, frame = 1.)
psicarrepsiprimecarre <- function(x)
{
basepsiprime(x) * basepsiprime(x) * basepsi(x) * basepsi(x) *
dnorm(x)
}
assign("psicarrepsiprimecarre", psicarrepsiprimecarre, frame = 1.)
psicarresurxcarre <- function(x)
{
(basepsi(x) * basepsi(x))/x^2. * dnorm(x)
}
assign("psicarresurxcarre", psicarresurxcarre, frame = 1.)
psiquatsurxcarre <- function(x)
{
(basepsi(x) * basepsi(x) * basepsi(x) * basepsi(x))/x^2. *
dnorm(x)
}
assign("psiquatsurxcarre", psiquatsurxcarre, frame = 1.)
epsiprime <- integrate(psiprime, lower = - Inf, upper = Inf)$integral
epsiprimecarre <- integrate(psiprimecarre, lower = - Inf, upper = Inf)$
integral
epsicarre <- integrate(psicarre, lower = - Inf, upper = Inf)$integral
epsiprimepsicarre <- integrate(psiprimepsicarre, lower = - Inf, upper
= Inf)$integral
epsicarrepsiprimecarre <- integrate(psicarrepsiprimecarre, lower = -
Inf, upper = Inf)$integral
epsicarresurxcarre <- integrate(psicarresurxcarre, lower = - Inf,
upper = Inf)$integral
epsiquatsurxcarre <- integrate(psiquatsurxcarre, lower = - Inf, upper
= Inf)$integral
return(epsiprime, epsiprimecarre, epsicarre, epsiprimepsicarre,
epsicarrepsiprimecarre, epsicarresurxcarre, epsiquatsurxcarre)
}
frob <- function(xx, yy, lambda, chuber = 1.345, n.orig = length(xx), compteurmax = 500., mytol = .Machine$double.eps^0.25, Smat = F)
{
# Fit of an M-type smoothing spline, as described in Eva Cantoni, Elvezio
# Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting Splines",
# Cahier du Département d'Econométrie, University of Geneva, 98.06, August 1998
if(length(xx) != length(yy)) stop("lengths of x and y must match")
# Computation of the weights and of the vector of unique values of "x" with the
# respective "y" values (cf. smooth.spline())
# The program works then with the sorted data (with respect to x)
if(missing(lambda)) stop("smoothing parameter lambda not supplied")
ww <- rep(1., length(xx))
nb <- length(xx)
sx <- unique(sort(xx))
myo <- match(xx, sx)
mynx <- length(sx)
if(nb != mynx) {
mycollaps <- .Fortran("suff",
as.integer(nb),
as.integer(mynx),
as.integer(myo),
as.double(xx),
as.double(yy),
as.double(ww),
double(mynx),
ybar = double(mynx),
wbar = double(mynx),
double(nb))
sy <- mycollaps$ybar
ww <- mycollaps$wbar
}
else {
sy <- yy[order(xx)]
}
esp.psi <- e.psi(chuber)
epsiprime <- esp.psi$epsiprime
epsicarre <- esp.psi$epsicarre
epsicarresurxcarre <- esp.psi$epsicarresurxcarre
epsiprimepsicarre <- esp.psi$epsiprimepsicarre
epsiquatsurxcarre <- esp.psi$epsiquatsurxcarre
epsiprimecarre <- esp.psi$epsiprimecarre
# Initialisation
epsicarrepsiprimecarre <- esp.psi$epsicarrepsiprimecarre
lin.fit <- lsfit(sx, sy)
frob.old <- sy - resid(lin.fit)
sigmahat.old <- scale.tau(frob.old, weights = ww)
myresid <- 1./sigmahat.old * resid(lin.fit)
compteur <- 0.
# Iterative solution for an M-type smoothing spline
coeff.spar <- (max(sx) - min(sx))^3.
# Cross-validation score. We suppose psi(0)=0.
repeat {
vecpsi <- psihuber(myresid, chuber)
ajust.new <- my.smooth.spline(sx, frob.old + sigmahat.old/
epsiprime * vecpsi, w = ww, n.orig = n.orig, spar = (
lambda * sigmahat.old)/epsiprime/coeff.spar/mean(ww),
all.knots = F)
frob.new <- ajust.new$y
diag.matsol <- ajust.new$lev
sigmahat.new <- scale.tau(sy, center = frob.new, weights = ww,
init.scale = sigmahat.old)
if(abs(max(frob.old - frob.new))/abs(max(myresid)) < mytol)
break
myresid <- 1./sigmahat.new * (sy - frob.new)
sigmahat.old <- sigmahat.new
frob.old <- frob.new
compteur <- compteur + 1.
if(compteur > compteurmax)
stop(paste("Convergence not obtained after",
compteurmax, "iterations"))
}
cv.score <- sigmahat.new^2./epsiprime^2. * weighted.mean((psihuber(
(yy - frob.new[myo])/sigmahat.new, chuber)/(1. - (diag.matsol[
myo] * rep(1., length(yy)))/ww[myo]))^2., rep(1., length(yy)))
sigma.ext <- NULL
RCp <- NULL
matsol <- NULL
# Output
if(Smat) {
Ident <- diag(1., nrow = length(sx))
# External estimation of scale for the RCp criterion
matsol <- matS(sx, (lambda * sigmahat.old)/epsiprime, wei = ww,
n.orig)
acoeff <- diff(sx)[-1.]/diff(sx, lag = 2.)
bcoeff <- diff(sx)[ - (length(sx) - 1.)]/diff(sx, lag = 2.)
normcoeff <- sqrt(ww[ - c(length(sy) - 1., length(sy))]^2. *
acoeff^2. + ww[ - c(1., 2.)]^2. * bcoeff^2. + ww[ -
c(1., length(sy))]^2.)
# RCp criterion
sigma.ext <- scale.tau((acoeff * ww[ - c(length(sy) - 1.,
length(sy))] * sy[ - c(length(sy) - 1., length(sy))] +
bcoeff * ww[ - c(1., 2.)] * sy[ - c(1., 2.)] - ww[
- c(1., length(sy))] * sy[ - c(1., length(sy))])/
normcoeff, center = 0.)
RCp.RASR <- sum(psihuber((yy - frob.new[myo])/sigma.ext, chuber
)^2.)
UmoinsV <- nb * epsicarre - 2./epsiprime * epsiprimepsicarre *
sum(diag(matsol)) + 1./sigma.ext^2. * (epsiprimecarre -
epsicarresurxcarre) * sum(diag((matsol - Ident) %*%
frob.new %*% t(frob.new) %*% t(matsol - Ident))) + 1./
epsiprime^2. * (epsicarrepsiprimecarre -
epsiquatsurxcarre) * sum(diag(matsol %*% t(matsol)))
RCp <- RCp.RASR - UmoinsV
}
estim <- frob.new
sigmahat <- sigmahat.new
x <- sx
yin <- sy
weights <- ww
Smatrix <- matsol
diagS <- diag.matsol
invisible(return(estim, sigmahat, sigma.ext, x, yin, myo, weights,
Smatrix, diagS, epsiprime, epsicarre, cv.score, RCp,
chuber, lambda))
}
matS <- function(xx, lambda, wei, n.orig)
{
# Computation of the matrix defining a resistant spline.
nb <- length(xx)
matS <- matrix(ncol = nb, nrow = nb)
Ident <- diag(1., nrow = nb)
myspar <- lambda/(max(xx) - min(xx))^3./mean(wei)
for(j in 1.:nb) {
matS[, j] <- my.smooth.spline(xx, Ident[, j], w = wei, n.orig
= n.orig, spar = myspar, all.knots = F)$y
}
matS
}
my.smooth.spline <- function(x, y, w, n.orig, df = 5., spar = 0., cv = F, all.knots = F, df.offset = 0., penalty = 1.)
{
#This function assumes that ties are treated before.
my.call <- match.call()
if(missing(y)) {
if(is.list(x)) {
if(any(is.na(match(c("x", "y"), names(x)))))
stop("cannot find x and y in list")
y <- x$y
x <- x$x
}
else if(is.complex(x)) {
y <- Im(x)
x <- Re(x)
}
else if(is.matrix(x) && ncol(x) == 2.) {
y <- x[, 2.]
x <- x[, 1.]
}
else {
y <- x
x <- time(x)
}
}
n <- length(x)
if(n != length(y))
stop("lengths of x and y must match")
if(missing(w))
w <- rep(1., n)
else {
if(length(x) != length(w))
stop("lengths of x and w must match")
#w <- (w * sum(w > 0))/sum(w) #this makes the weights sum to n
# Eva: This line was replaced by the following one, in order to construct
# the matrix S correctly, when there are ties in the x's (with matS.s)
if(any(w < 0.)) stop("all weights should be non-negative")
if(sum(w) != n.orig)
stop("sum of weigths muste be equal to n.orig")
}
if(missing(spar))
ispar <- 0.
else {
if(spar < 1.01e-15)
ispar <- 0.
else ispar <- 1.
}
if(cv) {
icrit <- 2.
}
else {
icrit <- 1.
}
dfinfo <- df.offset
if(!missing(df)) {
if(df > 1. & df < n) {
icrit <- 3.
dfinfo <- df
}
else warning("you must supply 1 < df < n")
}
n <- as.integer(n)
isetup <- as.integer(0.)
ier <- as.integer(1.)
mode(x) <- "double"
mode(y) <- "double"
mode(w) <- "double"
mode(spar) <- "double"
mode(ispar) <- "integer"
mode(icrit) <- "integer"
mode(isetup) <- "integer"
mode(ier) <- "integer"
x <- signif(x, 6.)
sx <- unique(sort(x))
o <- match(x, sx)
# Eva: this part is already done in frob.s
# collaps <- .Fortran("suff",
# n,
# as.integer(nx),
# as.integer(o),
# x,
# y,
# w,
# double(nx),
# ybar = double(nx),
# wbar = double(nx),
# double(n))
nx <- length(sx)
xbar <- (sx - sx[1.])/(sx[nx] - sx[1.])
if(all.knots) {
knot <- c(rep(xbar[1.], 3.), xbar, rep(xbar[nx], 3.))
nk <- nx + 2.
}
else {
knot <- .Fortran("sknotl",
xbar,
as.integer(nx),
knot = double(n + 6.),
k = integer(1.))
nk <- knot$k - 4.
knot <- knot$knot[seq(knot$k)]
}
low.parm <- 0.
high.parm <- 1.5
fit <- .Fortran("qsbart",
as.double(penalty),
as.double(dfinfo),
x = xbar,
y = y,
w = # y = collaps$ybar
w = w,
# w = collaps$wbar
as.integer(nx),
knot,
as.integer(nk),
coef = double(nk),
ty = double(nx),
lev = double(nx),
crit = double(1.),
c(icrit, ispar),
spar = spar,
c(as.double(low.parm), as.double(high.parm), as.double(0.001)),
isetup,
double((17. + nk) * nk),
as.integer(4.),
as.integer(1.),
ier = ier)
if(fit$ier > 0.) {
warning("smoothing parameter value too small or too large")
fit$ty <- rep(mean(y), nx)
}
lev <- fit$lev
df <- sum(lev)
# win <- collaps$wbar[1:nx]
win <- w[1.:nx]
if(cv) {
ww <- win
# for the next line
ww[!(ww > 0.)] <- 1.
cv.crit <- weighted.mean(((y - fit$ty[o])/(1. - (lev[o] * w)/
ww[o]))^2., w)
}
else cv.crit <- weighted.mean((y - fit$ty[o])^2., w)/(1. - (df.offset +
penalty * df)/sum(win))^2.
pen.crit <- sum(win * (y - fit$ty) * y)
fit.object <- list(knot = knot, nk = nk, min = sx[1.], range = sx[
nx] - sx[1.], coef = fit$coef)
oldClass(fit.object) <- "smooth.spline.fit"
object <- list(x = sx, y = fit$ty, w = win, yin = y[1.:nx], lev = #yin = collaps$ybar[1:nx]
lev = lev, cv.crit = cv.crit, pen.crit = pen.crit, df = df, spar = fit$
spar, fit = fit.object, call = my.call)
oldClass(object) <- "smooth.spline"
object
}
opt.RCp <- function(xx, yy, bornes.opt, ind.chuber = 1.345)
{
# Fit of an M-type smoothing spline with parameter chosen by robust C_p,
# as described in Eva Cantoni, Elvezio
# Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting Splines",
# Cahier du Département d'Econométrie, University of Geneva, 98.06, August 1998
if(length(xx) != length(yy)) stop("lengths of x and y must match")
# Expectation terms for the RCp criterion
nb <- length(xx)
psi.integr <- e.psi(ind.chuber)
epsiprime <- psi.integr$epsiprime
epsicarre <- psi.integr$epsicarre
epsicarresurxcarre <- psi.integr$epsicarresurxcarre
epsiprimepsicarre <- psi.integr$epsiprimepsicarre
epsiquatsurxcarre <- psi.integr$epsiquatsurxcarre
epsiprimecarre <- psi.integr$epsiprimecarre
epsicarrepsiprimecarre <- psi.integr$epsicarrepsiprimecarre
assign("ind.chuber", ind.chuber, frame = 1.)
assign("bornes.opt", bornes.opt, frame = 1.)
sx <- unique(sort(xx))
myo <- match(xx, sx)
Ident <- diag(1., length(sx))
# Objective function to minimize.
assign("Ident", Ident, frame = 1.)
# Minimisation
RCp.value <- function(lambda, xx, yy, nb, epsiprime, epsicarre,
epsiprimepsicarre, epsiprimecarre, epsicarresurxcarre,
epsicarrepsiprimecarre, epsiquatsurxcarre, myo)
{
ajust <- frob(xx, yy, lambda, chuber = ind.chuber, Smat = T)
fchap <- ajust$estim
Smatrix <- ajust$Smatrix
sigma.ext <- ajust$sigma.ext
ww <- ajust$weights
RCp.RASR <- sum(psihuber((yy - fchap[myo])/sigma.ext,
ind.chuber)^2.)
UmoinsV <- nb * epsicarre - 2./epsiprime * epsiprimepsicarre *
sum(diag(Smatrix)) + 1./sigma.ext^2. * (epsiprimecarre -
epsicarresurxcarre) * sum(diag((Smatrix - Ident) %*%
fchap %*% t(fchap) %*% t(Smatrix - Ident))) + 1./
epsiprime^2. * (epsicarrepsiprimecarre -
epsiquatsurxcarre) * sum(diag(Smatrix %*% t(Smatrix)))
RCp.RASR - UmoinsV
}
# Optimal values.
opt.lambda <- optimize(RCp.value, bornes.opt, xx = xx, yy = yy, nb = nb,
epsiprime = epsiprime, epsicarre = epsicarre, epsiprimepsicarre
= epsiprimepsicarre, epsiprimecarre = epsiprimecarre,
epsicarresurxcarre = epsicarresurxcarre, epsicarrepsiprimecarre
= epsicarrepsiprimecarre, epsiquatsurxcarre =
epsiquatsurxcarre, myo = myo)
lambda.opt <- opt.lambda$minimum
# Checks
RCp.opt <- opt.lambda$objective
mytol <- .Machine$double.eps^0.25
if(bornes.opt[2.] - lambda.opt < mytol)
warning(paste("lambda.opt=", round(lambda.opt, 6.),
"too close to upper bound", bornes.opt[2.],
"\n Search interval need to be widen"))
if(ind.chuber > 10000. && lambda.opt < mytol)
stop("lambda.opt cannot be obtained due to the instability of Cp. \n See comment of Figure 3, p. 145, in Cantoni & Ronchetti (2001)"
)
# Fit for the optimal value of the smoothing parameter
if(lambda.opt - bornes.opt[1.] < mytol) warning(paste("lambda.opt=",
round(lambda.opt, 6.), "too close to lower bound",
bornes.opt[1.], "\n Search interval need to be widen"))
# Output
ajust.opt <- frob(xx, yy, lambda.opt, chuber = ind.chuber, Smat = T)
invisible(return(lambda.opt, RCp.opt, ajust.opt))
}
opt.cv <- function(xx, yy, bornes.opt = c(0., 1.), ind.chuber = 1.345)
{
# Fit of an M-type smoothing spline with parameter chosen by robust
# cross-validation, as described in Eva Cantoni, Elvezio
# Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting Splines",
# Cahier du Département d'Econométrie, University of Geneva, 98.06, August 1998
if(length(xx) != length(yy)) stop("lengths of x and y must match")
nb <- length(xx)
psi.integr <- e.psi(ind.chuber)
epsiprime <- psi.integr$epsiprime
assign("ind.chuber", ind.chuber, frame = 1.)
# Distinct values of xx.
assign("bornes.opt", bornes.opt, frame = 1.)
sx <- unique(sort(xx))
# Objective function to minimize.
myo <- match(xx, sx)
# Minimisation
cv.value <- function(lambda, xx, yy, nb, epsiprime, myo)
{
ajust <- frob(xx, yy, lambda, chuber = ind.chuber)
fchap <- ajust$estim
diagS <- ajust$diagS
sigmahat <- ajust$sigmahat
ww <- ajust$weights
sigmahat^2./epsiprime^2. * weighted.mean((psihuber((yy - fchap[
myo])/sigmahat, ind.chuber)/(1. - (diagS[myo] * rep(
1., length(yy)))/ww[myo]))^2., rep(1., length(yy)))
}
# Optimal values.
opt.lambda <- optimize(cv.value, bornes.opt, xx = xx, yy = yy, nb = nb,
epsiprime = epsiprime, myo = myo)
lambda.opt <- opt.lambda$minimum
# Checks
cv.opt <- opt.lambda$objective
mytol <- .Machine$double.eps^0.25
if(bornes.opt[2.] - lambda.opt < mytol)
warning(paste("lambda.opt=", round(lambda.opt, 6.),
"too close to upper bound", bornes.opt[2.],
"\n Search interval need to be widen"))
# Fit for the optimal value of the smoothing parameter
if(lambda.opt - bornes.opt[1.] < mytol) warning(paste("lambda.opt=",
round(lambda.opt, 6.), "too close to lower bound",
bornes.opt[1.], "\n Search interval need to be widen"))
# Output
ajust.opt <- frob(xx, yy, lambda.opt, chuber = ind.chuber)
invisible(return(lambda.opt, cv.opt, ajust.opt))
}
psihuber <- function(x, chuber)
{
# This function defines the psi Huber function,
# with tuning constant equal to chuber.
x * pmin(1., chuber/abs(x))
}
Mspline/votfraud.dat 100644 20207 20154 450 7614007752 14113 0 ustar cantoni metri Machine Absentee
26427 346
15904 282
42448 223
19444 593
71797 572
-1017 -229
63406 671
15671 293
36276 360
36710 306
21848 401
65862 378
-13194 -829
56100 394
700 151
11529 -349
26047 160
44425 1329
45512 368
-5700 -434
51206 391
-564 1025
Mspline/simulated.example.dump 100644 20207 20154 11704 7514513561 16142 0 ustar cantoni metri "simex"<-
structure(.Data = list(xx = c(0.2159364684484899, 0.94926069444045424,
0.70554016530513763, 0.90520709287375212, 0.083439765963703394,
0.52705387957394123, 0.20753714488819242, 0.0012282119132578373,
0.84774063620716333, 0.5780115956440568, 0.56519387941807508,
0.82548357872292399, 0.30162268737331033, 0.270262086763978,
0.52639509830623865, 0.26024596486240625, 0.62108139554038644,
0.59766427427530289, 0.69887840142473578, 0.77176495688036084,
0.68722116248682141, 0.035150566603988409, 0.023296335246413946,
0.14422154519706964, 0.93533012829720974, 0.14017751393839717,
0.6998950382694602, 0.76077586458995938, 0.24018044583499432,
0.4084882540628314, 0.86981449043378234, 0.56475617224350572,
0.83972920989617705, 0.24705857690423727, 0.49881547596305609,
0.83144698711112142, 0.32386581506580114, 0.061407836619764566,
0.025746928527951241, 0.27463880646973848, 0.19253076333552599,
0.010002839379012585, 0.49294582707807422, 0.37892589950934052,
0.89308096375316381, 0.57220462337136269, 0.18636902514845133,
0.045602424535900354, 0.46904214564710855, 0.0056425123475492001,
0.067926475312560797, 0.198250079061836, 0.05044789332896471,
0.019665001425892115, 0.14864075183868408, 0.43060992751270533,
0.107445253059268, 0.84185505891218781, 0.62135014217346907,
0.96026840899139643, 0.045804535504430532, 0.70405057445168495,
0.43655805522575974, 0.6758022322319448, 0.85908479802310467,
0.34843307174742222, 0.28669627150520682, 0.98234925698488951,
0.23328682733699679, 0.8350418577902019, 0.46050783200189471,
0.96786748524755239, 0.30478319618850946, 0.92108155880123377,
0.16291052289307117, 0.96662822319194674, 0.53910678531974554,
0.41670937975868583, 0.88888543657958508, 0.34253728948533535,
0.51494773989543319, 0.041770535986870527, 0.96727905096486211,
0.028560718521475792, 0.30979224946349859, 0.2200382286682725,
0.31888188095763326, 0.80964374961331487, 0.62861587479710579,
0.84932251973077655, 0.39312346000224352, 0.38619101513177156,
0.90453245770186186, 0.25820326339453459, 0.88657401502132416,
0.66741975164040923, 0.60728075820952654, 0.074334654491394758,
0.13559750188142061, 0.072956295218318701), yy = c(4.3852782018970311,
-0.49613067470554584
, 0.51504749915218939, -0.58987956272749553, -1.2989375717931047,
0.77904270997032277, -0.55876262584276759, 0.0072072947863742205,
-0.34067126902621625, 0.72461791085987426, 1.2328083820594609,
-0.39746165940540112, -0.12080405014595533, -0.29122285483691451,
0.98452543329255726, -1.1728007869345913, 1.9010128910546955,
0.79806945561523535, 0.93778378063954004, 0.1994434589594333,
0.13417500559906564, -0.10166718099467176, -0.2896122137663274,
-1.0026986810985061, 0.46706143732706251, -1.4109768096475852,
0.65354831798715307, -0.36372184432238897, 5.5501757170904007,
1.2477516564608233, 0.75025885946248128, 0.9967340694660296,
0.85832157696453693, -0.72238360866601492, 0.14721069149627231,
-0.04188059576648176, 0.51976501853547141, -1.568596543175929,
-0.21809503611006703, -0.029494252976940805, -1.1041671157792587,
-0.35201097206679938, 0.91940756046474481, 0.2085002503044604,
0.0086350471268948625, 1.1708793119243095, -0.2854125542981053,
0.026526873417379249, 0.91280555641408856, -0.67319261219193782,
-0.83374123456251514, -0.10692042620272357, -0.90483262048909952,
0.063594640864058499, -0.5206663623040797, 0.82423707138925151,
-0.84850203645598576, -0.18756586796505847, 0.80112936129207191,
0.16900809154268628, 0.42784165458192513, -0.17459376326115872,
-0.34729437152242282, 0.66947033837440451, 0.463830116025779,
0.31360451196218053, 0.090195512157019736, -0.17589810138248582,
-0.83306794086052449, 0.31119926529776665, 0.24512753702790346,
-0.35726660030142654, 0.3319788832629636, -0.17321756414763273,
-0.82434455192664324, -0.060613917727520734, 1.3842070379510667,
1.8389362202431032, 0.67704902601748751, 0.43380332646341713,
0.69832416709017708, 0.068083092371948095, 0.70976044040617081,
-1.6926775284728146, 1.0670373992611364, -0.67186906135358404,
0.12973018031035607, 0.95513075225238786, 1.2648420881909352,
0.36276001946449177, 1.5389678693983353, -0.079245887189198339,
0.94070953776142463, 0.012478850240667583, -0.14733385176014058,
0.49864493356158573, 1.2347544783219293, -0.25859161101738637,
-1.0573573440473316, 2.4415726961017055)), row.names = c("1", "2", "3",
"4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28",
"29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40",
"41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52",
"53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64",
"65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76",
"77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88",
"89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100"
), class = "data.frame")
Mspline/help3.4/ 40755 20207 20154 0 7515017726 12672 5 ustar cantoni metri Mspline/help3.4/e.psi.d 100644 20207 20154 1601 7515017726 14150 0 ustar cantoni metri .BG
.FN e.psi
.TL
Auxiliary function for frob().
.DN
Auxiliary function for frob(), which computes the expectations terms (at the
normal model) appearing in the formulation of an M-type smoothing spline and in
the resistant criterion (RCV and RCp) for the selection of the smoothing parameter
.CS
e.psi(chuber=1.345)
.RA
.AG chuber
the tuning constant of the Huber psi function.
.RT
The function returns the expectations (at the normal model). See Eva Cantoni and
Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater for Smooting
Splines", 2001, Statistics and Computing, 11, 141-146, for more details.
.DT
The integrations are performed numerically via integrate().
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
frob, opt.cv, opt.RCp, integrate.
.WR
Mspline/help3.4/frob.d 100644 20207 20154 5552 7513557710 14073 0 ustar cantoni metri .BG
.FN frob
.TL
Fit an M-type smoothing spline to the input data.
.DN
Fit an M-type smoothing spline as described in Eva Cantoni and Elvezio Ronchetti,
"Resistant Selection of the Smoothing Paramater for Smooting Splines", 2001,
Statistics and Computing, 11, 141-146.
.CS
frob(xx, yy, lambda, chuber=1.345, n.orig=length(xx), compteurmax=500, mytol=.Machine$double.eps^0.25, Smat=F)
.RA
.AG xx
Values of the predictor variable.
.AG yy
Response variable, of the same length as xx.
.AG lambda
Smoothing parameter corresponding to the formulation of the reference paper.
Note that it is not the same parameter as spar in the function smooth.spline()
in S-PLUS. In fact spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
.OA
.AG chuber
Tuning constant of the Huber psi function. The default value is set to 1.345.
.AG n.orig
Number of original observations. This parameter controls the weighting procedure
when there are ties and is used to assess the condition sum(w)=n.orig in
my.smooth.spline(). The default is set to length(xx).
.AG compteurmax
Maximum number of iterations of the fitting algorithm.
.AG mytol
Convergence threshold.
.AG Smat
Logical value indicating wheter the whole smoother matrix has to be computed.
The default is F. In this case, only the diagonal elements of the smoother matrix
are computed. Moreover, the RCp criterion is not computed, and only the
cross-validation score is given.
.RT
A list is returned with the following components:
.AG estim
The fitted M-type smoothing spline corresponding to the ordered xx.
.AG sigmahat
The associated estimation of scale (Huber Proposal 2).
.AG sigma.ext
The external estimation of scale used in the construction of RCp (if Smat=T).
.AG x
The ordered distinct values of xx.
.AG yin
The y-values used at the unique x values (weighted averages of the input yy).
.AG myo
Value of match(xx, unique(sort(xx))). This allow one to recover the vector of
fitted values with respect to the original xx by considering estim[myo].
.AG weights
weights used in the fit. This has the same length as xx, and in the case of ties,
will consist of the accumulated weights at each unique value of x.
.AG Smatrix
The smoother matrix S, if Smat=T.
.AG diagS
The diagonal elements of the smoother matrix S.
.AG epsiprime
The expectation of the derivative of the psi function.
.AG epsicarre
The expectation of the squared psi function.
.AG cv.score
The value of the robust cross-validation criterion.
.AG RCp
The value of the robust Cp criterion, if Smat=T.
.AG chuber
The tuning constant used for the fit.
.AG lambda
The smoothing parameter used for the fit.
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
opt.cv, opt.RCp.
.EX
myfit<- frob(myx,myy,lambda=0.01)
motif()
plot(myx,myy)
lines(myfit$x,myfit$estim)
.KW robust, smooth.
.WR
Mspline/help3.4/matS.d 100644 20207 20154 1005 7513557720 14035 0 ustar cantoni metri .BG
.FN matS
.TL
Auxiliary function for frob().
.DN
Auxiliary function for frob() which computes the smoother matrix S.
.CS
matS(xx, lambda, wei, n.orig)
.RA
.AG xx
Values of the predictor variable.
.AG lambda
Smoothing parameter.
.AG wei
weights used in the fit.
.AG n.orig
Number of original observations.
.RT
The smoother matrix S.
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
frob.
.WR
Mspline/help3.4/my.smooth.spline.d 100644 20207 20154 2113 7513557752 16365 0 ustar cantoni metri .BG
.FN my.smooth.spline
.TL
Auxiliary function for frob().
.DN
Auxiliary function for frob(). It is a modified version of smooth.spline, which
assumes that ties are treated before. (This is done in the first part of frob()).
.CS
my.smooth.spline(x, y, w, n.orig, df=5, spar=0, cv=F, all.knots=F, df.offset=0,
penalty=1)
.RA
.AG x
Values of the predictor variable.
.AG y
Response variable, of the same length as x.
.AG w
Weights used in the fit.
.AG n.orig
Number of original observations.
.AG spar
Smoothing parameter satisfying spar = (lambda * sigma)/epsiprime/(max(xx) - min(xx))^3.
.OA
The arguments of my.smooth.spline listed below are in general not used. They are
present because my.smooth.spline is a modified version of smooth.spline.
.AG df
.AG cv
.AG all.knots
.AG df.offset
.AG penalty
.RT
This function returns an object of class smooth.spline. See the help of this
function for further details.
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
frob.
.WR
Mspline/help3.4/opt.RCp.d 100644 20207 20154 2627 7513557776 14444 0 ustar cantoni metri .BG
.FN opt.RCp
.TL
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust Cp.
.DN
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
Cp as described in Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of
the Smoothing Paramater for Smooting Splines", 2001, Statistics and Computing,
11, 141-146.
.CS
opt.RCp(xx, yy, bornes.opt, ind.chuber=1.345)
.RA
.AG xx
Values of the predictor variable.
.AG yy
Response variable, of the same length as xx.
.AG bornes.opt
Interval in which the optimal smoothing parameter has to be found.
.OA
.AG ind.chuber
Tuning constant of the Huber psi function. The default value is set to 1.345.
.RT
A list is returned with the following components:
.AG lambda.opt
The optimal value of the smoothing parameter obtained by minimization of the
robust Cp criterion.
.AG cv.opt
The value of the robust Cp for the optimal value of the smoothing parameter.
.AG ajust.opt
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
frob,opt.cv.
.EX
myoptimalfit<- opt.RCp(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
.KW robust, smooth.
.WR
Mspline/help3.4/opt.cv.d 100644 20207 20154 2715 7514766723 14361 0 ustar cantoni metri .BG
.FN opt.cv
.TL
Fit an M-type smoothing spline to the input data with the smoothing parameter
chosen by robust cross-validation.
.DN
Fit an M-type smoothing spline with the smoothing parameter chosen by robust
cross-validation as described in Eva Cantoni and Elvezio Ronchetti, "Resistant
Selection of the Smoothing Paramater for Smooting Splines", 2001, Statistics and
Computing, 11, 141-146.
.CS
opt.cv(xx, yy, bornes.opt, ind.chuber=1.345)
.RA
.AG xx
Values of the predictor variable.
.AG yy
Response variable, of the same length as xx.
.AG bornes.opt
Interval in which the optimal smoothing parameter has to be found.
.OA
.AG ind.chuber
Tuning constant of the Huber psi function. The default value is set to 1.345.
.RT
A list is returned with the following components:
.AG lambda.opt
The optimal value of the smoothing parameter obtained by minimization of the
robust cross-validation criterion.
.AG cv.opt
The value of the robust cross-validation for the optimal value of the smoothing
parameter.
.AG ajust.opt
The fit for the optimal value of the smoothing parameter. It is an object of type
frob(). See this function for more details.
.SH REFERENCES
Eva Cantoni and Elvezio Ronchetti, "Resistant Selection of the Smoothing Paramater
for Smooting Splines", 2001, Statistics and Computing, 11, 141-146.
.SA
frob,opt.RCp.
.EX
myoptimalfit<- opt.cv(myx,myy,c(0,1))
motif()
plot(myx,myy)
lines(myoptimalfit$ajust.opt$x,myoptimalfit$ajust.opt$estim)
.KW robust, smooth.
.WR
Mspline/help3.4/psihuber.d 100644 20207 20154 362 7513560052 14727 0 ustar cantoni metri .BG
.FN psihuber
.TL
Huber psi function.
.DN
Define a Huber psi function.
.CS
psihuber(x, chuber)
.RA
.AG x
Value at which the function needs to be evaluated.
.AG chuber
Value of the tuning constant.
.RT
Value of the function at point x.
.WR