Przewidywanie odpowiedzi z nowych krzywych za pomocą pakietu fda w języku R.

Zasadniczo wszystko, co chcę zrobić, to przewidzieć odpowiedź skalarną za pomocą niektórych krzywych. Dotarłem do regresji (używając fRegress z pakietu fda), ale nie mam pojęcia, jak zastosować wyniki do NOWEGO zestawu krzywych (do przewidywania).

Mam N = 536 krzywych i 536 odpowiedzi skalarnych. Oto co dotychczas zrobiłem:

• Stworzyłem podstawę dla krzywych.
• Stworzyłem obiekt fdPar w celu wprowadzenia kary
• Stworzyłem obiekt fd za pomocą smooth.basis, aby wygładzić krzywe wybraną karą na określonej podstawie.
• Przeprowadziłem regresję za pomocą fRegress (), regresując krzywe odpowiedzi skalarnej.

Teraz wszystko, co chciałbym zrobić, to użyć tej regresji do stworzenia prognoz dla nowego zestawu danych, które mam. Nie mogę znaleźć łatwego sposobu na zrobienie tego.

Twoje zdrowie

Nie dbam o fda„s wykorzystania Incepcja -Jak list-wewnątrz-listy-w całym wykazie struktur obiektów, ale moja odpowiedź będzie trwać przez system autorzy pakietów zostały utworzone.

Myślę, że warto najpierw pomyśleć o tym, co dokładnie robimy. W oparciu o opis tego, co zrobiłeś do tej pory, robię to, co według mnie robisz (daj mi znać, jeśli coś źle zinterpretowałem). Będę nadal używać notacji, a ze względu na brak rzeczywistych danych, przykład z analizy danych funkcjonalnych Ramsaya i Silvermana oraz analizy danych funkcjonalnych Ramsaya, Hookera i Gravesa za pomocą R i MATLAB (niektóre z poniższych równań i kodu są podnoszone bezpośrednio z tych książek).

Modelujemy odpowiedź skalarną za pomocą funkcjonalnego modelu liniowego, tj

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

W pewnym stopniu rozszerzamy . Używamy, powiedzmy, funkcji bazowych K. Więc,$\beta$$K$

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

W notacji macierzowej jest to .$\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$

Rozszerzamy również funkcje współzmienne w pewnym stopniu (powiedzmy funkcje podstawowe ). Więc,$L$

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

Ponownie, w notacji macierzowej jest to .$X(s)=\mathbf{C}\boldsymbol{\psi}(s)$

A zatem, jeśli pozwolimy , nasz model można wyrazić jako$\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$

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

$\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$$\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$

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

I wygląda nam to znacznie lepiej.

Teraz widzę, że dodajesz jakąś regularyzację. fdaPakiet współpracuje z chropowatości kar formularza

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

$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)$

$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}$

dlatego naszym problemem jest jedynie regresja kalenicy z rozwiązaniem:

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

Przeszedłem powyższe, ponieważ: (1) myślę, że ważne jest, abyśmy rozumieli, co robimy, i (2) niektóre z powyższych są konieczne, aby zrozumieć część kodu, którego użyję później. Do kodu ...

Oto przykład danych z kodem R. Korzystam z kanadyjskiego zestawu danych pogodowych zawartych w fdapakiecie. Modelujemy dzienne opady atmosferyczne dla wielu stacji pogodowych za pomocą funkcjonalnego modelu liniowego i wykorzystujemy profile temperatur (temperatury rejestrowano raz dziennie przez 365 dni) z każdej stacji jako funkcjonalne zmienne towarzyszące. Postępujemy podobnie do tego, co opisujesz w swojej sytuacji. Dane rejestrowano na 35 stacjach. Podzielę zestaw danych na 34 stacje, które zostaną wykorzystane jako moje dane, oraz na ostatnią stację, która będzie moim „nowym” zestawem danych.

Kontynuuję przez kod R i komentarze (zakładam, że znasz się na fdapakiecie tak, że nic poniżej nie jest zbyt zaskakujące - jeśli tak nie jest, daj mi znać):

# 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)

Teraz, gdy rok temu dowiedziałem się o danych funkcjonalnych, bawiłem się tym pakietem. Nie mogłem też predict.fRegressdać mi tego, czego chciałem. Patrząc na to teraz, wciąż nie wiem, jak to zrobić. Będziemy musieli więc uzyskać prognozy półautomatycznie. Będę używał elementów, dla których wyciągnąłem prosto z kodu fRegress(). Ponownie kontynuuję za pomocą kodu i komentarzy.

Po pierwsze, konfiguracja:

# 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)

Teraz, aby uzyskać prognozy

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

fRegressyhatfdobjfRegressyhatfdobj$\int_0^T X_i (s) \beta(s)$$X_i$$\beta$

Zwykle fRegressoblicza dopasowane wartości, zapętlając zmienne towarzyszące zapisane w annPrecTemp$xfdlist. Dlatego w przypadku naszego problemu zastępujemy tę listę zmiennych towarzyszących odpowiadającą jej na naszej nowej liście zmiennych, tj templistNew. Oto kod (identyczny z kodem znalezionym w fRegressdwóch edycjach, usunięciu niepotrzebnego kodu i dodaniu kilku komentarzy):

# 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
}

(uwaga: jeśli spojrzysz na ten fragment i otaczający go kod fRegress, zobaczysz kroki, które przedstawiłem powyżej).

Przetestowałem kod, ponownie uruchamiając przykład pogody, używając wszystkich 35 stacji jako naszych danych i porównałem wynik z powyższej pętli do annPrecTemp$yhatfdobji wszystko pasuje. Uruchomiłem go również kilka razy, używając różnych stacji jako moich „nowych” danych i wszystko wydaje się rozsądne.

Daj mi znać, jeśli którekolwiek z powyższych jest niejasne lub jeśli coś nie działa poprawnie. Przepraszamy za zbyt szczegółową odpowiedź. Nie mogłem się powstrzymać :) A jeśli jeszcze ich nie masz, sprawdź dwie książki, w których pisałem tę odpowiedź. To naprawdę dobre książki.