Ich verwende Paket fda
insbesondere Funktion fRegress
. Diese Funktion enthält eine weitere Funktion, die eigchk
heißt und prüft, ob die Koeffizientenmatrix singulär ist.Erfassen der Druck der Funktion
Hier ist die Funktion, wie die Paketbesitzer (J. O. Ramsay, Giles Hooker und Spencer Graves) es geschrieben haben.
eigchk <- function(Cmat) {
# check Cmat for singularity
eigval <- eigen(Cmat)$values
ncoef <- length(eigval)
if (eigval[ncoef] < 0) {
neig <- min(length(eigval),10)
cat("\nSmallest eigenvalues:\n")
print(eigval[(ncoef-neig+1):ncoef])
cat("\nLargest eigenvalues:\n")
print(eigval[1:neig])
stop("Negative eigenvalue of coefficient matrix.")
}
if (eigval[ncoef] == 0) stop("Zero eigenvalue of coefficient matrix.")
logcondition <- log10(eigval[1]) - log10(eigval[ncoef])
if (logcondition > 12) {
warning("Near singularity in coefficient matrix.")
cat(paste("\nLog10 Eigenvalues range from\n",
log10(eigval[ncoef])," to ",log10(eigval[1]),"\n"))
}
}
Als Sie das letzte Mal sehen, wenn die Bedingung überprüft, ob logcondition
größer als 12 und druckt dann die Bereiche der Eigenwerte.
Der folgende Code implementiert die Verwendung der Regularisierung mit Rauheitspennalty. Der Code stammt aus dem Buch "Funktionelle Datenanalyse mit R und Matlab".
annualprec = log10(apply(daily$precav,2,sum))
tempbasis =create.fourier.basis(c(0,365),65)
tempSmooth=smooth.basis(day.5,daily$tempav,tempbasis)
tempfd =tempSmooth$fd
templist = vector("list",2)
templist[[1]] = rep(1,35)
templist[[2]] = tempfd
conbasis = create.constant.basis(c(0,365))
betalist = vector("list",2)
betalist[[1]] = conbasis
SSE = sum((annualprec - mean(annualprec))^2)
Lcoef = c(0,(2*pi/365)^2,0)
harmaccelLfd = vec2Lfd(Lcoef, c(0,365))
betabasis = create.fourier.basis(c(0, 365), 35)
lambda = 10^12.5
betafdPar = fdPar(betabasis, harmaccelLfd, lambda)
betalist[[2]] = betafdPar
annPrecTemp = fRegress(annualprec, templist, betalist)
betaestlist2 = annPrecTemp$betaestlist
annualprechat2 = annPrecTemp$yhatfdobj
SSE1.2 = sum((annualprec-annualprechat2)^2)
RSQ2 = (SSE - SSE1.2)/SSE
Fratio2 = ((SSE-SSE1.2)/3.7)/(SSE1/30.3)
resid = annualprec - annualprechat2
SigmaE. = sum(resid^2)/(35-annPrecTemp$df)
SigmaE = SigmaE.*diag(rep(1,35))
y2cMap = tempSmooth$y2cMap
stderrList = fRegress.stderr(annPrecTemp, y2cMap, SigmaE)
betafdPar = betaestlist2[[2]]
betafd = betafdPar$fd
betastderrList = stderrList$betastderrlist
betastderrfd = betastderrList[[2]]
Als Strafe Faktor verwenden die Autoren bestimmte lambda
. Der folgende Code implementiert die Suche nach dem passenden `Lambda.
loglam = seq(5,15,0.5)
nlam = length(loglam)
SSE.CV = matrix(0,nlam,1)
for (ilam in 1:nlam) {
lambda = 10ˆloglam[ilam]
betalisti = betalist
betafdPar2 = betalisti[[2]]
betafdPar2$lambda = lambda
betalisti[[2]] = betafdPar2
fRegi = fRegress.CV(annualprec, templist,
betalisti)
SSE.CV[ilam] = fRegi$SSE.CV
}
Durch den Wert des loglam
und Kreuzvalidierung zu ändern Ich nehme die besten lambda
equaire, doch wenn die Länge des loglam
zu groß ist oder seine Werte der Koeffizientenmatrix zu singulrity führen. Ich erhalte die folgende Meldung:
Log10 Eigenvalues range from
-5.44495317739048 to 6.78194912518214
durch die Funktion Erstellt eigchk
wie ich bereits oben erwähnt habe. Jetzt ist meine Frage, gibt es eine Möglichkeit, diese so genannte Warnung zu fangen? Mit Fang meine ich eine Funktion oder Methode, die mich warnt, wenn dies passiert ist und ich die Werte der loglam
anpassen konnte. Da es in der Funktion neben diesem Ausdruck der Nachricht keine eigentliche Warndefinition gibt, ging mir die Idee aus.
Vielen Dank für Ihre Anregungen.
Wenn das eine Nachricht ist, können Sie es mit 'tryCatch' fangen, aber ich bin nicht sicher, was Sie danach tun möchten – rawr
@raw Mein Ziel ist eine Funktion, die die Werte des Vektors 'Loglam' basiert auf diesem Druck. Aber das Problem ist, dass dieser Druck eine Art zusätzlicher Hinweis und keine echte Warnung ist. – user3185925