2016-04-05 22 views
0

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.

+0

Wenn das eine Nachricht ist, können Sie es mit 'tryCatch' fangen, aber ich bin nicht sicher, was Sie danach tun möchten – rawr

+0

@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

Antwort

0

Durch "Fang die Warnung", wenn Sie meinen, wird Sie darauf aufmerksam machen, dass es ein mögliches Problem mit loglam gibt, dann möchten Sie vielleicht try und tryCatch Funktionen betrachten. Dann können Sie das Verhalten definieren, das Sie implementieren möchten, wenn eine Warnbedingung erfüllt ist.

Wenn Sie nur die Ausgabe der Warnung speichern möchten (was vom Titel der Frage angenommen werden könnte, aber möglicherweise nicht das ist, was Sie wollen), dann schauen Sie sich capture.output an.

+0

Ich habe mir die tryCatch-Funktion bereits angesehen und wenn mein Verständnis stimmt, wäre dies der Fall, wenn meine Funktion die Warnung als Ausgabe hätte, aber dieser Ausdruck ist eine Art zusätzlicher Hinweis, der Evaluierungen nicht stoppt. – user3185925

Verwandte Themen