From a course by Prof. Camelia Goga
April 24, 2024
library(dda)
library(tidyverse)
library(ggplot2)
data(vietnam_temperature)
vntemp <- vietnam_temperature |>
filter(year == 2016) |>
mutate(t_max = lapply(t_max, unclass), id = seq_len(n())) |>
unnest_longer(t_max, indices_to = "day")
vntemp |> ggplot(aes(day, t_max, group = id, color = region)) +
geom_point() +
geom_line() +
xlab("Day number") +
ylab("Temperature (°C)")
Figure 1: There are more discretisation points (365 or 366) than observations (63 provinces).
Figure 2: The discretisation points are not equally spaced and different for each curve. Source: CONTENT study, Sixth Framework Programme of the European Union.
Figure 3: The x-axis is not time.
We observe a sample of random functions Y1,…,Yn: Yi:Ω⟶F and we consider that Yi,i=1,…,n are iid.
We always observe the functions Yi,i=1,…,n in a finite (often large) number of discretisation points (which may be different from one curve to another): ti1,…,tiDi⊂[0,T] so the data is presented in vector form.
To overcome these difficulties,
In functional statistics, the objects of interest are functions: Y=(Y(t))t∈[0,T]
Therefore, in theory, Y could be measured at any time t∈[0,T].
Nevertheless, in practice:
To overcome these difficulties, the true curve Y=(Y(t))t∈[0,T] can be “reconstructed” from the discretised data using non-parametric smoothing methods.
library(fda)
y <- filter(vietnam_temperature, year == 2016 & province == "HANOI") |>
pull(t_max) |>
unlist()
breaks <- seq(1, length(y), length.out = 10)
bsp <- create.bspline.basis(breaks = breaks)
yhat <- smooth.basis(seq_len(length(y)), y, bsp)
invisible(
plot(yhat, xlab = "Day number", ylab = "Temperature (°C)")
)
points(seq_len(length(y)), y)Figure 4: Smoothing of daily maximum temperature in 2016 in Hanoi.
We can use non-parametric methods like
to give an estimate ˆf(td) of f(t);
The re-constructed curve is Y(t)=ˆf(t) for all t∈[0,T].
Let b1(⋅),…,bq(⋅) be a basis of spline functions (B-splines or truncated polynomials) of dimension q=K+m with K equidistant interior knots and piecewise polynomials of degree m−1;
We denote by bT(t)=(b1(t),…,bq(t)) the function vectors of the basis calculated at a point t∈[0,T] and by B the matrix containing the values of bT(t) for the discretisation points td,d=1,…,D: B=(bT(t1)bT(t2)…bT(tD))=(b1(t1)…bq(t1)b1(t2)…bq(t2)………b1(tD)…bq(tD))
The estimator ˆf(td) of f(td) is given by \begin{aligned} \widehat f(t_d) &= \mathbf{b}^T(t_d)\boldsymbol{\widehat{\theta}} \quad d=1, \ldots, D\quad \mbox{with}\\ \boldsymbol{\widehat{\theta}} &= \underset{\boldsymbol{\theta}\in \mathbb R^q}{\operatorname{argmin}} \sum_{d=1}^D(y_d-\mathbf{b}^T(t_d)\boldsymbol{\theta})^2=(\mathbf{B}^T\mathbf{B})^{-1}\mathbf{B}^T\mathbf{Y}_d \end{aligned} where \mathbf{Y}_d=(Y_d)_{d=1}^D.
In the end, we have the curve reconstructed by non-parametric regression using spline functions: \widehat Y(t)=\widehat f(t)=\mathbf{b}^T(t)\boldsymbol{\widehat\theta}=\mathbf{a}^T(t)\mathbf{Y}_d \quad \text{for all}\quad t\in [0, T];
The number of knots K can be chosen using the cross-validation or generalised cross-validation criterion and m is small (3 or 4). The knots can be positioned at the quantiles of the instants t_1, \dots, t_D.
We have a sample of n curves, Y_1, \ldots, Y_n which are, for simplicity, measured at the same discretisation times t_1, \ldots, t_D;
Let Y_i(t_d)=Y_{id} for i=1, \ldots, n and d=1, \ldots, D.
Let n be non-parametric models: \begin{aligned} Y_{1d} & = f_1(t_d)+\varepsilon_{1d}, \quad d=1, \ldots, D \\ Y_{2d} & = f_2(t_d)+\varepsilon_{2d}, \quad d=1, \ldots, D \\ &\vdots \\ Y_{nd} & = f_n(t_d)+\varepsilon_{nd}, \quad d=1, \ldots, D \end{aligned}
The objective is to estimate f_i, i=1, \ldots, n;
Let b_1(\cdot), \ldots, b_q(\cdot) be the basis of spline functions of order q=K+m: K equidistant interior knots of degree m-1;
This same basis is used to re-construct all the n curves;
Using the method described for smoothing a single curve, we obtain: \begin{aligned} \widehat Y_i(t)= \widehat f_i(t) &= \mathbf{b}^T(t)\widehat{\boldsymbol{\theta}}_i \\ &= \mathbf{b}^T(t)(\mathbf{B}^T\mathbf{B})^{-1}\mathbf{B}^T\mathbf{Y}_i, \quad i=1, \ldots, n \end{aligned} where \mathbf{Y}_i=(Y_{id})_{d=1}^D;
The same method can be used in the case where the discretisation times are not the same for all the curves. In this situation, we consider the union of all the measurement instants and the knots are placed at the quantiles of this new set.
The variables Y_1, \ldots, Y_n are now functions of t\in [0, T]: Y_i=(Y_i(t))_{t\in [0, T]}, and we consider that they are independent and identically distributed;
We consider that the Y_i, i=1, \ldots, n belong to the Hilbert space L^2[0, T], the space of integrable square functions with the usual scalar product and the associated norm: \langle f,g \rangle = \int_0^Tf(t)g(t)dt, \quad \|f \| =\sqrt{\int_0^Tf^2(t) \mathrm dt}
The processes (functions) Y_i are defined on a probability space (\Omega, \mathcal K, \mathbb P) with values in L^2[0, T]: Y_i(\omega)=(Y_i(t,\omega))_{t\in [0, T]}, \quad i=1, \ldots, n
Remark. This is not the Lebesgue integral but rather the Bochner integral (its generalisation to Banach spaces).
Let Y_i, i=1, \ldots, n be independent and identically distributed with Y.
Figure 5: Empirical mean of the NIR curves, computed pointwise.
Under the preceding conditions: \mathbb E(\langle Y, f \rangle)=\langle \mathbb E(Y), f\rangle \quad \mbox{for all}\quad f\in L^2[0, T]; so the integral and the scalar product commute.
“Pythagorean theorem”: \mathbb E\|Y-\mu\|^2=\mathbb E\|Y\|^2-\mu^2
The \Gamma operator has the following properties (assuming that \mathbb E(\|Y\|^2)<\infty):
\displaystyle \langle \Gamma f, g\rangle=\mathbb E\left(\langle Y-\mu, f\rangle\langle Y-\mu, g\rangle\right) for all f, g\in L^2[0, T];
\Gamma=\mathbb E (Y\otimes Y)-\mu\otimes\mu
it is a self-adjoint operator (\Gamma=\Gamma^*, where \Gamma^* is the adjoint operator) because the kernel \gamma(\cdot, \cdot) is symmetric and non-negative (\langle \Gamma u, u \rangle \geq 0);
it has finite trace: \|\Gamma\|_{TR}=\sum_{j\geq 1}\langle (\Gamma^*\Gamma)^{1/2} e_j, e_j\rangle=\mathbb E\|Y\|^2<\infty;
is a Hilbert-Schmidt operator: |\Gamma\|^2_{HS}=\sum_{j\geq 1}\|\Gamma e_j\|^2<\infty,
where \{e_j\}_{j\geq 1} is any orthonormal basis of L^2[0, T].
Let (\lambda_j, v_j)_{j\leq 1} with \lambda_j\in \mathbb R and v_j=(v_j(t))_{t\in [0, T]}\in L^2[0, T] be called eigenvalues and associated eigenvectors (functions) of \Gamma: \begin{aligned} \Gamma v_j &= \lambda_j \ v_j\\ \Gamma v_j (t)& = \lambda_j \ v_j(t), \quad t \in [0,T], \quad j\geq 1 \end{aligned}
\Gamma is non-negative so \lambda_j\geq 0 for all j\geq 1;
\Gamma is self-adjoint, so the eigenvectors associated with the distinct eigenvalues are linearly independent and orthogonal;
The operator \Gamma is compact and self-adjoint so the set of non-zero eigenvalues is finite or is a sequence which converges to zero. Let (\lambda_j)_{j\geq 1} be the ordered eigenvalues: \lambda_1\geq\lambda_2\geq\ldots>0. The associated eigenvectors (v_j)_{j\geq 1} form an orthonormal basis in L^2[0,T] for \overline{\operatorname{Im} (\gamma)} i.e \langle v_j, v_{j'} \rangle = \delta_{jj'}. The covariance then has the following decomposition: \Gamma=\sum_{j\geq 1}\lambda_jv_j\otimes v_j \quad \mbox{et}\quad \Gamma f=\sum_{j\geq 1}\lambda_j \langle f,v_j\rangle v_j;
The covariance function \gamma associated with \Gamma has the following decomposition: \gamma(s,t)=\sum_{j\geq 1}\lambda_jv_j(s) v_j(t) \quad \mbox{for any}\quad s,t\in [0, T];
Y-\mu=\sum_{j\geq 1}\langle Y-\mu, v_j\rangle v_j and \mathbb E(\langle Y-\mu, v_j\rangle v_j)=0;
The trace norm of \Gamma is given by \|\Gamma\|_{TR}=\sum_{j\geq 1}\langle \Gamma v_j, v_j\rangle=\sum_{j\geq 1} \lambda_j<\infty and the Hilbert-Schmidt norm by \|\Gamma\|_{HS}=\sum_{j\geq 1} \lambda^2_j<\infty.
\Gamma^{1/2}=\sum_{j\geq 1}\lambda_j^{1/2}v_j\otimes v_j and \Gamma^p=\sum_{j\geq 1}\lambda_j^pv_j\otimes v_j.
Therefore, the optimal subspace of dimension q, (unique if \lambda_j is distinct), is the space generated by the q eigenvectors associated to the q largest eigenvalues of \Gamma;%Thus the optimal subspace of dimension q, which is unique if \lambda_q > \lambda_{q+1}, is the space generated by the q eigenfunctions of \Gamma associated to the q largest eigenvalues. Y_k(t) = \mu(t) + \sum_{j =1}^q \langle Y_k - \mu,v_j \rangle v_j(t) + R_{q,k}(t), \quad t \in [0,T]
The space generated by v_1, \cdots, v_q gives a representation of the variation over time of the data around the mean.
The variance of the projection on v_j is given by: \lambda_j = \frac{1}{n} \sum_{k=1}^n \langle Y_k - \mu,v_j \rangle^2 \ because \frac{1}{n} \sum_{k=1}^n \langle Y_k - \mu,v_j \rangle=0.
The notion of quantile is not generalisable to \mathbb R^p because we cannot give a total order relation in \mathbb R^p; the same thing in a more general space like L^2[0, T];
Several definitions have been proposed of a concept that could be considered a quantile but none has the basic property of a quantile of order \alpha in R, having \alpha\% of the distribution less than this value;
A definition for the median m = \underset{y\in L^2[0,T]}{\operatorname{argmin}} \sum _{i=1}^n \|Y_i - y\|;
Let X_1, \ldots, X_p be real variables and \mathbf{x}'_i=(1,X_{i1}, \ldots, X_{ip});
The data are (\mathbf{x}'_i, Y_i(t)) for i=1, \ldots, n;
The following linear model (Faraway 1997) M1:\quad Y_i(t)=\mathbf{x}'_i\boldsymbol{\beta}(t)+\varepsilon_{i}(t), \quad t\in [0,T] where \boldsymbol{\beta}(t)=(\beta_0(t), \beta_1(t), \ldots, \beta_p(t))' is the functional regression coefficient and the functional residuals \varepsilon_{k}(t) are centred : \mathbb E(\varepsilon_{i}(t))=0 and decorrelated with each the same covariance function \gamma: \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{j}(t))=0, \quad i\neq j, t\in [0,T] \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{i}(r))=\gamma(t,r), \quad t,r\in [0,T]
It is assumed that the matrix \mathbf{X}=(\mathbf{x}'_i)_{i=1}^n of size (n,p+1) is of rank p+1;
The pointwise estimator of \boldsymbol{\beta} (for t fixed) is obtained by ordinary least squares: \begin{aligned} \widehat{\boldsymbol{\beta}}(t) & = \underset{\beta(t)}{\operatorname{argmin}} \sum_{i=1}^n(Y_i(t)-\mathbf{x}'_i\boldsymbol{\beta}(t))^2=\|\mathbf{Y}(t)-\mathbf{X}\boldsymbol{\beta}(t)\|^2 \\ & = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}\mathbf{Y}(t) \end{aligned} where \mathbf{Y}(t)=(Y_i(t))_{i=1}^n;
It can be shown that the estimator \widehat{\boldsymbol{\beta}}=(\widehat{\boldsymbol{\beta}}(t))_{t\in [0, T]} satisfies the global criterion: \widehat{\boldsymbol{\beta}}=\mbox{argmin}_{\beta}\int_{0}^T\|\mathbf{Y}(t)-\mathbf{X}\boldsymbol{\beta}(t)\|^2dt
The predictions are \widehat{\mathbf{Y}}(t)=\mathbf{X}\widehat{\boldsymbol{\beta}}(t)=\mathbf{P}_X\mathbf{Y}(t);
The estimated residuals are \widehat{\boldsymbol{\varepsilon}}(t)=\mathbf{Y}(t)-\widehat{\mathbf{Y}}(t)=\left(\mathbf{I}_n-\mathbf{P}_X\right)\mathbf{Y}(t): \mathbb E(\widehat{\boldsymbol{\varepsilon}}(t))=0 \operatorname{cov}(\widehat{\boldsymbol{\varepsilon}}(t), \widehat{\boldsymbol{\varepsilon}}(r))=\gamma(t,r)\left(\mathbf{I}_n-\mathbf{P}_X\right) \mathbb E(\widehat{\boldsymbol{\varepsilon}}(t)'\widehat{\boldsymbol{\varepsilon}}(r))=(n-p-1)\gamma(t,r) so an unbiased estimator of \gamma(t,r) is given by \widehat{\gamma}(t,r)=\frac{1}{n-p-1}\sum_{i=1}^n(Y_i(t)-\widehat Y_i(t))(Y_i(r)-\widehat Y_i(r))
R
Figure 7: Implementation of scalar-on-function regression with refund::pfr.
Let X_1, \ldots, X_p be functional variables, X_j=(X_j(t))_{t\in [0, T]}; let \mathbf{x}'_i(t)=(1,X_{i1}(t), \ldots, X_{ip}(t));
The data are (\mathbf{x}'_i(t), Y_i(t)) for i=1, \ldots, n; M1:\quad Y_i(t)=\mathbf{x}'_i(t)\boldsymbol{\beta}(t)+\varepsilon_{i}(t), \quad t\in [0,T] where \boldsymbol{\beta}(t)=(\beta_0(t), \beta_1(t), \ldots, \beta_p(t))' is the functional regression coefficient and the functional residuals \varepsilon_{k}(t) are centred: \mathbb E(\varepsilon_{i}(t))=0 and decorrelated with each the same covariance function \gamma: \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{j}(t))=0, \quad i\neq j, t\in [0,T] \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{i}(r))=\gamma(t,r), \quad t,r\in [0,T]
It is assumed that the matrix \mathbf{X}(t)=(\mathbf{x}'_i(t))_{i=1}^n of size (n,p+1) is of rank p+1 for each t\in [0, T];
The point estimator of \boldsymbol{\beta} (for t fixed) is obtained by ordinary least squares: \begin{aligned} \widehat{\boldsymbol{\beta}}(t) & = \underset{\beta(t)}{\operatorname{argmin}} \sum_{i=1}^n(Y_i(t)-\mathbf{x}'_i(t)\boldsymbol{\beta}(t))^2=\|\mathbf{Y}(t)-\mathbf{X}(t)\boldsymbol{\beta}(t)\|^2\\ & = (\mathbf{X(t)}'\mathbf{X}(t))^{-1}\mathbf{X}(t)\mathbf{Y}(t) \end{aligned} where \mathbf{Y}(t)=(Y_i(t))_{i=1}^n;
It can be shown that the estimator \widehat{\boldsymbol{\beta}}=(\widehat{\boldsymbol{\beta}}(t))_{t\in [0, T]} satisfies the global criterion: \widehat{\boldsymbol{\beta}}= \underset{\beta}{\operatorname{argmin}} \int_{0}^T\|\mathbf{Y}(t)-\mathbf{X}(t)\boldsymbol{\beta}(t)\|^2 \mathrm dt
