Vorhersage der Reaktion neuer Kurven mit dem FDA-Paket in R.

Grundsätzlich möchte ich nur anhand einiger Kurven eine skalare Reaktion vorhersagen. Ich bin so weit gekommen, eine Regression durchzuführen (mit fRegress aus dem FDA-Paket), habe aber keine Ahnung, wie ich die Ergebnisse auf einen NEUEN Satz von Kurven anwenden soll (zur Vorhersage).

Ich habe N = 536 Kurven und 536 skalare Antworten. Folgendes habe ich bisher getan:

• Ich habe eine Basis für die Kurven erstellt.
• Ich habe ein fdPar-Objekt erstellt, um eine Strafe einzuführen
• Ich habe das fd-Objekt mit glatt.basis erstellt, um die Kurven mit der gewählten Strafe auf der angegebenen Basis zu glätten.
• Ich habe mit fRegress () eine Regression ausgeführt und dabei die Kurven der Skalarantwort zurückgeführt.

Jetzt möchte ich nur noch diese Regression verwenden, um Vorhersagen für einen neuen Datensatz zu erstellen, den ich habe. Ich kann anscheinend keinen einfachen Weg finden, dies zu tun.

Prost

Die fdaVerwendung von Inception- ähnlichen Objektstrukturen für Listen innerhalb von Listen innerhalb von Listen ist mir egal , aber meine Antwort wird sich an das System halten, das die Paketschreiber erstellt haben.

Ich finde es lehrreich, zuerst darüber nachzudenken, was wir genau tun. Basierend auf Ihrer Beschreibung dessen, was Sie bisher getan haben, glaube ich, dass Sie dies tun (lassen Sie mich wissen, wenn ich etwas falsch interpretiert habe). Ich werde weiterhin die Notation verwenden und aufgrund fehlender realer Daten ein Beispiel aus der Funktionsdatenanalyse von Ramsay und Silverman sowie die Funktionsdatenanalyse von Ramsay, Hooker und Graves mit R und MATLAB (Einige der folgenden Gleichungen und Codes werden direkt aufgehoben aus diesen Büchern).

Wir modellieren eine skalare Antwort über ein funktionales lineares Modell, dh

$y_i = \beta_0 + \int_0^T X_i(s)\beta(s)ds + \epsilon_i$

Wir erweitern die in gewisser Weise. Wir verwenden beispielsweise Basisfunktionen. Damit,$\beta$$K$

$\beta(s) = \sum\limits_{k=1}^K b_k\theta_k (s)$

In der Matrixnotation ist dies .$\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$

Wir erweitern auch die kovariaten Funktionen in einigen Basen (sagen wir Basisfunktionen). Damit,$L$

$X_i(s) = \sum\limits_{k=1}^L c_{ik}\psi_k (s)$

In der Matrixnotation ist dies wiederum .$X(s)=\mathbf{C}\boldsymbol{\psi}(s)$

Wenn wir also , kann unser Modell ausgedrückt werden als$\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$

$y = \beta_0 + \mathbf{CJb}$ .

Und wenn wir und , ist unser Modell$\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$$\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$

$\mathbf{y} = \mathbf{Z}\boldsymbol{\xi}$

Und das kommt uns viel vertrauter vor.

Jetzt sehe ich, dass Sie eine Art Regularisierung hinzufügen. Das fdaPaket arbeitet mit Rauheitsstrafen des Formulars

$P = \lambda\int[L\beta(s)]^2 ds$

für einige linearen Differentialoperator . Jetzt kann gezeigt werden (Details hier weggelassen - es ist wirklich nicht schwer, dies zu zeigen), dass, wenn wir die Strafmatrix als definieren$L$$\mathbf{R}$

$\mathbf{R}= \lambda\left(\begin{array}{cccc} 0 & 0 & \cdots & 0 \\ 0 & R_1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_K \end{array}\right)$

wo in Bezug auf die Basis Ausbau der ist , dann minimieren wir die bestraft Summe der Quadrate:$R_i$$\beta_i$

$\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right)'\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right) + \lambda\boldsymbol{\xi}'\mathbf{R}\boldsymbol{\xi}$ ,

und so ist unser Problem nur eine Gratregression mit Lösung:

$\hat{\boldsymbol{\xi}}=(\mathbf{Z}'\mathbf{Z}+\lambda\mathbf{R})^{-1}\mathbf{Z}'\mathbf{y}$ .

Ich habe das oben Gesagte durchgearbeitet, weil (1) ich denke, es ist wichtig, dass wir verstehen, was wir tun, und (2) einige der oben genannten Dinge notwendig sind, um einen Teil des Codes zu verstehen, den ich später verwenden werde. Weiter zum Code ...

Hier ist ein Datenbeispiel mit R-Code. Ich verwende den im fdaPaket enthaltenen kanadischen Wetterdatensatz . Wir modellieren den logarithmischen Jahresniederschlag für eine Reihe von Wetterstationen über ein funktionelles lineares Modell und verwenden Temperaturprofile (die Temperaturen wurden 365 Tage lang einmal täglich aufgezeichnet) von jeder Station als funktionale Kovariaten. Wir werden ähnlich vorgehen, wie Sie es in Ihrer Situation beschreiben. Die Daten wurden an 35 Stationen aufgezeichnet. Ich werde den Datensatz in 34 Stationen aufteilen, die als meine Daten verwendet werden, und die letzte Station, die mein "neuer" Datensatz sein wird.

Ich fahre mit R-Code und Kommentaren fort (ich gehe davon aus, dass Sie mit dem fdaPaket so vertraut sind, dass nichts im Folgenden zu überraschend ist - wenn dies nicht der Fall ist, lassen Sie es mich bitte wissen):

# pick out data and 'new data'
dailydat <- daily$precav[,2:35]
dailytemp <- daily$tempav[,2:35]
dailydatNew <- daily$precav[,1]
dailytempNew <- daily$tempav[,1]

# set up response variable
annualprec <- log10(apply(dailydat,2,sum))

# create basis objects for and smooth covariate functions
tempbasis <- create.fourier.basis(c(0,365),65)
tempSmooth <- smooth.basis(day.5,dailytemp,tempbasis)
tempfd <- tempSmooth$fd

# create design matrix object
templist <- vector("list",2)
templist[[1]] <- rep(1,34)
templist[[2]] <- tempfd

# create constant basis (for intercept) and
# fourier basis objects for remaining betas
conbasis <- create.constant.basis(c(0,365))
betabasis <- create.fourier.basis(c(0,365),35)
betalist <- vector("list",2)
betalist[[1]] <- conbasis
betalist[[2]] <- betabasis

# set roughness penalty for betas 
Lcoef <- c(0,(2*pi/365)^2,0)
harmaccelLfd <- vec2Lfd(Lcoef, c(0,365))
lambda <- 10^12.5
betafdPar <- fdPar(betabasis, harmaccelLfd, lambda)
betalist[[2]] <- betafdPar

# regress
annPrecTemp <- fRegress(annualprec, templist, betalist)

Als ich vor ungefähr einem Jahr zum ersten Mal über funktionale Daten unterrichtet wurde, habe ich mit diesem Paket herumgespielt. Ich konnte predict.fRegressmir auch nicht geben, was ich wollte. Wenn ich jetzt zurückblicke, weiß ich immer noch nicht, wie ich es dazu bringen soll, sich zu verhalten. Wir müssen die Vorhersagen also nur halb manuell abrufen. Ich werde Teile verwenden, für die ich direkt aus dem Code gezogen habe fRegress(). Wieder gehe ich über Code und Kommentare weiter.

Zunächst das Setup:

# create basis objects for and smooth covariate functions for new data
tempSmoothNew <- smooth.basis(day.5,dailytempNew,tempbasis)
tempfdNew <- tempSmoothNew$fd

# create design matrix object for new data
templistNew <- vector("list",2)
templistNew[[1]] <- rep(1,1)
templistNew[[2]] <- tempfdNew

# convert the intercept into an fd object
onebasis <- create.constant.basis(c(0,365))
templistNew[[1]] <- fd(matrix(templistNew[[1]],1,1), onebasis)

Nun zu den Vorhersagen

$\hat{\mathbf{y}}_{\mathrm{new}}=\mathbf{Z}_{\mathrm{new}}\hat{\boldsymbol{\xi}}$

Ich nehme nur den Code, fRegressmit dem er berechnet yhatfdobjund leicht bearbeitet wird. fRegressberechnet yhatfdobjdurch Schätzen des Integrals über die Trapezregel (wobei und in ihren jeweiligen Basen erweitert sind). $\int_0^T X_i (s) \beta(s)$$X_i$$\beta$

Normalerweise fRegresswerden die angepassten Werte berechnet, indem die in gespeicherten Kovariaten durchlaufen werden annPrecTemp$xfdlist. Für unser Problem ersetzen wir diese Kovariatenliste durch die entsprechende in unserer neuen Kovariatenliste, d templistNew. H. Hier ist der Code (identisch mit dem Code, der fRegressmit zwei Änderungen, einigen Löschungen von nicht benötigtem Code und einigen Kommentaren hinzugefügt wurde):

# set up yhat matrix (in our case it's 1x1)
yhatmat <- matrix(0,1,1)

# loop through covariates
p <- length(templistNew)
for(j in 1:p){
    xfdj       <- templistNew[[j]]
    xbasis     <- xfdj$basis
    xnbasis    <- xbasis$nbasis
    xrng       <- xbasis$rangeval
    nfine      <- max(501,10*xnbasis+1)
    tfine      <- seq(xrng[1], xrng[2], len=nfine)
    deltat     <- tfine[2]-tfine[1]
    xmat       <- eval.fd(tfine, xfdj)
    betafdParj <- annPrecTemp$betaestlist[[j]]
    betafdj    <- betafdParj$fd
    betamat    <- eval.fd(tfine, betafdj)
    # estimate int(x*beta) via trapezoid rule
    fitj       <- deltat*(crossprod(xmat,betamat) - 
                      0.5*(outer(xmat[1,],betamat[1,]) +
              outer(xmat[nfine,],betamat[nfine,])))
    yhatmat    <- yhatmat + fitj
}

(Hinweis: Wenn Sie sich diesen Block und den umgebenden Code in fRegressansehen, werden die oben beschriebenen Schritte angezeigt.)

Ich habe den Code getestet, indem ich das Wetterbeispiel mit allen 35 Stationen als Daten erneut ausgeführt und die Ausgabe der obigen Schleife mit verglichen habe annPrecTemp$yhatfdobjund alles übereinstimmt. Ich habe es auch ein paar Mal mit verschiedenen Stationen als "neue" Daten ausgeführt und alles scheint vernünftig.

Lassen Sie mich wissen, wenn einer der oben genannten Punkte unklar ist oder wenn etwas nicht richtig funktioniert. Entschuldigen Sie die zu detaillierte Antwort. Ich konnte mir nicht helfen :) Und wenn Sie sie noch nicht besitzen, lesen Sie die beiden Bücher, mit denen ich diese Antwort geschrieben habe. Das sind wirklich gute Bücher.