| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
Hilbert space embeddings of distributions, or kernel mean embeddings, map distributions into reproducing kernel Hilbert spaces (RKHS’s) where kernel methods are extended to probability measures. Kernel mean embedding methods have many applications in probabilistic modeling of complex phenomena, statistical inference based on high-dimensional data, and artificial intelligence. In this TCIU Reproducing Kernel Hilbert Spaces and Temporal Distribution Dynamics (TDD) Section we explore the relations between RKHS, kime-phase distributions, and temporal distribution dynamics.
Specifically, in this section we explore repeated spatiotemporal sampling, where \(\forall\ n\), we observe \(K\) (longitudinal or spatiotemporal) sample means \(\{\bar{X}_{n1},\bar{X}_{n2},\cdots, \bar{X}_{nK}\}\). In this case, \(\forall\ 1\leq k_o\leq K\) each IID sample \(\{X_{1k_o},X_{2k_o},\cdots ,X_{nk_o}\}\sim F_X\) (across the kime-phase space) allows us to compute the \(k_o^{th}\) sample mean \(\bar{X}_{nk_o}\).
Computing \(\bar{X}_{nk_o}\) from the \(k_o^{th}\) sample \(\{X_{1k_o},\cdots ,X_{nk_o}\}\) is identical to directly sampling \(\bar{X}_{nk_o}\) from the sampling distribution of the mean \(F_{\bar{X}_n}\). Kime-phase sampling, at a fixed spatiotemporal index \(k=k_o\) samples from the kime-phase distribution, \(\{X_{1k_o},\cdots , X_{nk_o}\}\). The complementary spatiotemporal sampling reflects the varying \(k\) index, where the kime-phase variability is annihilated by using a kime-phase aggregator, in this case sample-averages, to generate a sample \(\{\bar{X}_{n1},\cdots , \bar{X}_{nK} \}\).
Prior to data acquisition, \(\bar{X}_{nk_o}\) is a random variable, once the observed data values are plugged in, it’s a constant. Hence, the sample mean random variable (a sum of \(n\) \(\sim F_X\) variables) \(\bar{X}_{nk_o}\) based on \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim F_X\), and the (single!) random variable \(Y\sim F_{\bar{X}_{nk_o}}\) represent exactly the same random variable.
There is a kime-phase connection to quantum fluctuations, statistical inference, and artificial intelligence. This is related to random drawing of a second-order tensor \(X_{n\times K}\) (a random matrix) representing \(K\) (spatiotemporal or longitudinal) samples of \(n\) IID observations \(\{X_{1k_o},\cdots , X_{nk_o}\}_{k_o}\sim F_X\). Computing the \(K\) corresponding means \(\bar{X}_{nk_o}=\frac{1}{n}\sum_{i=1}^n{X_{ik_o}}\) is equivalent to randomly drawing \(K\) samples (a random vector) directly from the sampling distribution of the mean \(F_{\bar{X}_{nk_o}}\).
Kernel functions are positive definite operators that arise as a way to perform inner products \(\langle \cdot, \cdot \rangle\) in high-dimensional feature spaces \(\mathcal{F}\) given observed data points \(x, y \in \mathcal{X}\). These kernel functions guarantee the existence of mappings from the native state spaces of observed data \(\mathcal{X}\) to dot product feature spaces \(\mathcal{F}\), \(\phi: \mathcal{X} \to \mathcal{F}\), such that the kernel inner product \(\kappa(x, y) = \langle \phi(x), \phi(y)\rangle\) can be computed without the need for explicit calculation of the enigmatic feature mapping \(\phi\).
Kernel methods are used by machine learning (ML) algorithms and artificial intelligence (AI) techniques to replace the native space inner products \(\langle x, y\rangle_{\mathcal{X}}\) with other, possibly nonlinear, distance or similarity measures computed in the higher dimensional feature space representing the range (output) of nonlinear transformations \(\phi: \mathcal{X} \to \mathcal{F}\), which map the native space (domain, \(\mathcal{X}\)) into a higher-dimensional feature space \(\mathcal{F}\). Each feature map \(\phi\) transforms the native space inner product distance calculations \(\langle x, y\rangle_{\mathcal{X}}: \mathcal{X}\times \mathcal{X}\to \mathbb{R}\) with kernel function calculations \(\kappa(x, y) = \langle \phi(x), \phi(y)\rangle_{\mathcal{F}}: \mathcal{F}\times \mathcal{F}\to \mathbb{R}\) on the feature space \(\mathcal{F}\). As most of the time the elements of the inner product are clearly identifiable, we can skip the space-designating subscript and write \(\langle \cdot, \cdot\rangle\) for \(\langle \cdot, \cdot\rangle_{\mathcal{X}}\) and \(\langle \cdot, \cdot\rangle_{\mathcal{F}}\).
In other words, we interpret \(\kappa(x, y)\) as a nonlinear similarity measure between the sets \(x\) and \(y\). Many ML/AI algorithms that depend only on the inner product (distances) between state space data points, i.e., on \(\langle x, y\rangle\). Hence, this kernel-trick allows us to obtain nonlinear extensions of such algorithms by simply substituting the kernel distance measures \(\kappa(x, y) = \langle \phi(x), \phi(y)\rangle_{\mathcal{F}}\) for the corresponding native space distances, \(\langle x, y\rangle\) without altering the algorithm learning, training, and validation steps. For instance, consider a polynomial feature map \(\phi:\mathbb{R}^2\to \mathbb{R}^3\)
\[\phi(x) = (x^2_1,\ x^2_2,\ \sqrt{2}\ x_1 x_2),\ \forall\ x=(x_1,x_2)\in\mathbb{R}^2\ ,\]
which transforms the calculations of the native-space inner product to another feature-space inner product as follows
\[\underbrace{\langle x, y\rangle_{\mathcal{X}} \equiv x_1 y_1 + x_2 y_2}_{\text{native space similarity, }\mathcal{X}} \longrightarrow \underbrace{\kappa(x, y) = \langle \phi(x), \phi(y)\rangle_{\mathcal{F}}\equiv x_1^2 y_1^2 + x_2^2 y_2^2 + 2x_1 x_2 y_1 y_2}_{\text{(corresponding) RKHS space similarity, }\mathcal{F}}\equiv \langle x, y\rangle_{\mathcal{X}}^2 \ .\]
In this case, the kernel trick suggests that the \(\phi\)-transformed similarity measure, \(\langle \phi(x), \phi(y)\rangle_{\mathcal{F}}\), is just the square of the inner product in the native space \(\mathcal{X}\),i.e., \(\langle x, y\rangle ^2\).
The more general \(d\)-degree polynomial kernel relies on a map \(\phi(x):\mathcal{X}\equiv\mathbb{R}^n \to \mathcal{F}\), where \(\phi(x)\) represent all possible \(d^{th}\) degree ordered products of the observable data \(x,y\in\mathbb{R}^n\)
\[\kappa(x, y) = \langle \phi(x), \phi(y)\rangle_{\mathcal{F}} = \langle x, y\rangle_{\mathcal{X}} ^d\ .\]
In DSPA Chapter 6 we show several types of kernels commonly used in support vector machines (SVM) and other ML methods. Here are some kernel examples.
\[\kappa(x, y)=\langle x,\ y \rangle\ .\]
\[\kappa(x, y) = \langle x, y\rangle ^d=(x\cdot y+1)^d \ .\]
\[\kappa(x, y)=\tanh(s\ \langle x, y\rangle + o).\]
\[\kappa(x, y)=e^{- \left (\frac{\lVert x-y\rVert^2}{2\sigma^2}\right )} .\]
\[\kappa(x, y)=\tanh(s\langle x, y\rangle + o)\ .\]
\[\kappa(x, y) = e^{-\left (\frac{\lVert x-y\rVert}{\sigma}\right )} \ . \]
\[\kappa_{\nu} (x,y)=\sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)}\left (\sqrt{2\nu}\frac{\|x-y\|}{\rho}\right )^{\nu} K_{\nu}\left (\sqrt{2\nu}\frac{\|x-y\|}{\rho} \right )\ ,\]
where \(\sigma^2\) is the variance parameter and \(\rho\) controls the kernel scale (width).
When \(\nu=\frac{1}{2}\), \[\kappa_{\frac{1}{2}}(x,y)=\sigma^2\sqrt{\frac{2}{\pi}\frac{\|x-y\|}{\rho}} \underbrace{K_{\frac{1}{2}}\left (\frac{\|x-y\|}{\rho}\right )}_{\left (2\frac{\|x-y\|}{\rho\pi}\right )^{-1/2} \exp\left (-\frac{\|x-y\|}{\rho}\right )}= \sigma^2 \exp\left (-\frac{\|x-y\|}{\rho}\right )\ .\]
As \(\nu\to\infty\), the Matern kernel transforms into the Gaussian kernel,
\[\lim_{\nu\to\infty}\kappa_{\nu}(x,y)=\sigma^2\exp \left (-\frac{\|x-y\|^2}{\rho}\right )\]
where the Laplace method may be used to approximate the asymptotic convergence of the kernel as\(\nu\to\infty\)
\[K_{\nu}(x)\to \sqrt{\frac{\pi}{2\nu}}e^{-\nu}(2\nu/x)^{\nu},\ \Gamma(\nu)\to \sqrt{2\pi}\nu^{\nu-1/2}e^{-\nu}\ ,\]
\[K_{\nu}(\sqrt{v}x)\to \sqrt{\frac{\pi}{2\nu}}e^{-\nu}\left (2\frac{\sqrt{\nu}}{x}\right )^{\nu}\exp(x^2)\ ,\]
and \(\sqrt{\nu x^2+\nu^2}\sim \nu+\frac{x^2}{2},\ \nu\to \infty\).
This unique interpolation of the Matern kernel between Gaussian and Laplacian kernels is effectively a 1-parameter kernel homotopy, just like the t-distribution is a 1-parameter (probability distribution) homotopy between Cauchy distribution and Gaussian distribution.
Also, Matern kernels generate classical Sobolev spaces, indicating that reproducing-kernel Hilbert spaces (RKHS) associated with the Green function can be embedded into or are equivalent to a generalized Sobolev space.
Here are the Matern kernels for a couple of \(\nu\in \left \{\frac{3}{2}, \frac{5}{2} \right \}\) parameters:
\[\kappa_{\frac{3}{2}}(x,y)=\sigma^2\left (1+\frac{\sqrt{3}\|x-y\|}{\rho}\right ) \exp\left (-\frac{\sqrt{3}\|x-y\|}{\rho}\right)\]
\(\kappa_{\frac{5}{2}}(x,y)=\sigma^2(1+\frac{\sqrt{5}\|x-y\|}{\rho}+\frac{5(x-y)^2}{3\rho^2})\exp(-\frac{\sqrt{3}\|x-y\|}{\rho})\)
Below are plots of the shapes of several Matern kernels.
library(GPBayes)
besselk_nu <- function(nu) {
return(function(x) {
return(GPBayes::BesselK(nu, x))
})
}
maternkernel <- function(sigma,nu,rho,x){
scaled_x = sqrt(2*nu)*x/rho
return (sapply(scaled_x,besselk_nu(nu))*scaled_x**nu*sigma**2*2**(1-nu)/gamma(nu))
}
library(plotly)
p <- plot_ly(type="scatter", mode="lines")
x_values = seq(from = -5, to = 5, by = 0.01)
p <- p %>%add_trace(x=~x_values, y=~maternkernel(1,1/8,1,abs(x_values)),
name = ~paste("nu=", 1/4), mode="lines")
p <- p %>%add_trace( x = ~x_values, y = ~maternkernel(1,1/2,1,abs(x_values)),
name = ~paste("nu=", 1/2), mode="lines")
p <- p %>%add_trace( x = ~x_values, y = ~maternkernel(1,1,1,abs(x_values)),
name = ~paste("nu=", 1), mode="lines")
p <- p %>%add_trace( x = ~x_values, y = ~maternkernel(1,5,1,abs(x_values)),
name = ~paste("nu=", 5), mode="lines")
p <- p %>%add_trace(x = ~x_values, y = ~maternkernel(1,50,1,abs(x_values)),
name = ~paste("nu=", 50), mode="lines")
p = layout(p, title = "Matern kernels for different parameter values",
xaxis = list(title="X values"), yaxis = list(title="Y values"))
p\[\kappa(x, y) =−Bessel_{\nu +1}^{n} \left ( ∥x−y∥^2 /\sigma^2 \right ) \ ,\]
\[Bessel_{\nu +1}^{n}(u) = \frac{1}{(4 \pi)^{(\nu +1)/2}\Gamma ((\nu +1)/2)} \int_0^{\infty} {\frac{e^{-\left (\frac{\pi \vert u \vert^2}{s}+\frac{s}{4 \pi}\right )}} {s^{1 + \frac{n - (\nu + 1)}{2}}}\,\mathrm{d}s}\ , \forall\ u\in\mathbb{R}^n\setminus\{0\} .\]
\[\kappa(x, y) = \prod_{i=1}^n\left (1+ x_i y_i (1+\min(x_i,y_i))- \frac{(x_i+y_i)\min^2(x_i,y_i)}{2} +\frac{\min^3(x_i,y_i)}{3}\right )\ .\]
The R implementations of all kernlab kernels are available here.
In general, for each pair of observations \(x,y\in\mathcal{X}\), evaluating the kernel \(\kappa(x, y)\) involves two computationally expensive steps:
Note that there is a simplified alternative way to compute directly the kernel inner product (second step) \(\langle \phi(x), \phi(y)\rangle _{\mathcal{F}}\) without explicitly computing the feature maps (step one). This requires an explicit relation, like the (quadratic and \(d^{th}\) degree) polynomial kernel relations we saw earlier, e.g., \(\kappa(x, y) = \langle x, y\rangle_{\mathcal{F}}^2\) or \(\kappa(x, y) = \langle x, y\rangle_{\mathcal{F}}^d\).
It turns out that for each positive definite kernel \(\kappa(x, y)\), there always exists a feature map \(\phi: \mathcal{X} \to \mathcal{F}\) so \(\kappa(x, y) = \langle \phi(x), \phi(y)\rangle _{\mathcal{F}}\). And since feature space inner product \(\langle \cdot, \cdot\rangle _{\mathcal{F}}\) is positive definite, it implies that \(\kappa(x, y) = \langle \phi(x), \phi(y)\rangle _{\mathcal{F}}\) is positive definite \(\forall\ \phi:\mathcal{X} \to \mathcal{F}\).
Definition: A reproducing kernel Hilbert space (RKHS) \(\mathcal{H}\) is defined as the space of feature functions \(\phi:\mathcal{X} \to \mathcal{F}\), corresponding to reproducing positive semi-definite and symmetric kernels \(\kappa\)
\[\kappa: \mathcal{X}\times \mathcal{X} \longrightarrow \mathbb{C}\ \ ({\text{or }} \mathbb{R})\ ,\]
\[\text{Reproducing}: f(x)=\langle \kappa(x,\cdot),f(\cdot)\rangle_{\mathcal{H}} \forall f\in\mathcal{H}\] \[\text{Positive semi-definiteness}:\sum_{i,j=1}^n {c_i^*\ c_j\ \kappa(x_i,x_j)}\geq 0,\ \forall\ n\in\mathbb{N},\ \forall\ \{x_i\}_i \in \mathcal{X}, \forall\ \{c_i\}_i \in \mathbb{C} \ . \]
Two important RKHS properties are:
\[\langle \kappa (x, \cdot ),\ \kappa (x_o, \cdot )\rangle_{\mathcal{H}}\equiv \kappa (x, x_o)\ .\]
Thus, we do not need to explicate the feature map \(\phi_x = \kappa(x, \cdot )\) as it can be derived directly from the kernel \(\kappa\) by the above reproducing relation.
Note: Below we will demonstrate several examples of using kernel/RKHS approaches for regularized linear modeling of data (Ridge and LASSO regression) and show the notion of LASSO shrinkage on features (identification of salient variables associated with predicting a specific outcome of interest) and shrinkage on cases (selection of observation instances, i.e., participants, that play the most influential role in deriving model-based predictions of the outcome).
Let’s look at some hands-on parameter estimation examples of RKHS including a regularized linear model fitting (Ridge and LASSO) and simulation cases.
Let’s consider the DSPA amyotrophic lateral sclerosis (ALS) case-study, which consists of \(2,223\) observations and \(131\) numeric variables. We select ALSFRS slope as our outcome variable, as it captures the patients’ clinical decline over a year. In this simple demonstration, we’ll only use \(5\) covariates \(\{x_i\}_{i=1}^5\), \(x_1=\)Age_mean, \(x_2=\)Albumin_max, \(x_3=\)Albumin_median, \(x_4=\)Albumin_min, \(x_5=\)Albumin_range, to predict the outcome \(Y=ALSFRS\_slope\). This ALS case-study is intentionally chosen as it represents a tough modeling, prediction, and classification AI challenge.
We start by importing the raw data and identifying the design matrix \(X\) and the clinical outcome of interest, \(Y=ALSFRS\_slope\).
library(tidyverse)
# devtools::install_github("kupietz/kableExtra")
# See: https://stackoverflow.com/questions/76118194/error-when-loading-kableextra-in-markdown-file-after-updating-to-r-4-3-0
library(kableExtra)
ALS.train <- read.csv("https://umich.instructure.com/files/1789624/download?download_frd=1")
# summary(ALS.train); colnames(ALS.train)
y = ALS.train$ALSFRS_slope
X = as.matrix(ALS.train[ , 2:6])
X = apply(X, 2, scales::rescale, to = c(0, 1))To avoid potential computationally singular results, we define a new function inverse(), which will be used in place of the more commonly used base::solve().
Define a bivariate reproducing linear kernel operator, rk(),
\(\kappa(s,t) = rk(s,t)=s'\cdot t=\langle s, t \rangle\), in this case a simple inner-product.
For a set of vectors \(S = \{x_1,x_2,\cdots, x_n\}\) representing an observed sample of \(n\) multivariable data points \(\{x_i\}\), the Gram matrix is defined as the kernel-matrix
\[G_{n\times n}=\left (G_{ij} \equiv \langle \phi(x_i), \phi(x_j) \rangle \equiv \underbrace{\kappa (x_i, x_j )}_{kernel} \right )\ ,\]
where the kernel function \(\kappa(\cdot,\cdot)\) computes the inner products in the feature-space based on the feature map \(\phi(\cdot)\). Note that \(G_{ij}\equiv G_{ji},\ \forall\ 1\leq i,j\leq n\) and the symmetric Gram matrix, \(G'\equiv G\), contains all the information needed to compute the pairwise distances within the dataset. Indeed, the Gram matrix is a lossy representation of the complete original dataset (set of features). Despite the fact that the Gram matrix loses some information about the original dataset it has some desirable properties:
The power of the Gram matrix comes from its preservation of all pairwise distances between the observed data points, i.e., the vectors \(\{x_i\}_{i=1}^n\), and their \(\phi\)-transformed counterparts, \(\{\phi(x_i)\}_{i=1}^n\).
In our examples, the gram(X,RK) function will represent an alternative of the standard (square transposed cross-product matrix) \(tcrossprod(X)=tcrossprod(X,Y\equiv X)= X\times Y'\) function, which takes the (transposed) cross-product of a matrix with its transpose, i.e., \(X\times X'\). The gram() operator allows us to generalize tcrossprod(), so that it can be used with alternative kernels (RK) provided as arguments to gram().
Recall the two alternative forms of the generalized cross-validation (GCV) score defined in terms of the root sum-square (RSS) error or the deviance
\[GCV_{RSS} = \frac{RSS/n}{(1 − \nu/n)^2} \\ GCV_{Deviance} = -2\frac{LogLikelihood/n}{(1 − \nu/n)^2} \ ,\]
where \(\nu\) is the effective degree of freedom. Here we will define the ridge(X, y, lambda) function, which computes the Gram matrix \(G_{n\times n}\), formulates a non-singular design matrix \(Q\) (with an intercept first column),
computes the inverse of the cross-product matrix, \(M\), estimates the ridge parameters, and calculates the GCV measure.
The R package glmnet already provides an efficient method for fitting elastic net regularized linear models.
We will define a de novo ridge() method that demonstrates RKHS approaches using
the Gram matrix to fit a Ridge regularized linear model for estimating the parameters (effect sizes).
Recall that if \(X\) is the design matrix of observed values of predicting-covariates and \(Y\) is the observed outcome vector, the linear model predictions (fitted-values) are expressed as
\[\underbrace{\hat{Y}}_{fitted\ values}=X\hat{\beta}=X(X'X)^{-1}X'Y\ , \]
where the effect-size estimates \(\hat{\beta}=\{\beta_i\}_i\) for all covariates \(\{x_i\}_i\) correspond to the specific objective (cost/loss) function used in the optimization process
We implement the kernelized ridge regression as follows (Some additional details See appendix 5.3):
ridge <- function(X, y, lambda,kernel=rk) {
K = gram(X,kernel) #assumed X to be n*d where observations stacked on the rows, GramMatrix (nxn) XTX (already phi^Tphi)
rows = dim(X)[1]
n = length(y)
#Q = cbind(1, GramMatrix) no need for intercept as the kernel already captures this
#S = rbind(0, cbind(0, GramMatrix))
M = K + lambda*diag(dim(K)[1]) # (singularity protection)K+lambda I)^{-1}
M_inv = inverse(M) # invert M
omega_tilde = M_inv %*% y
f_hat = K %*% omega_tilde
A = K %*% M_inv # Hat matrix maps response to fitted values
tr_A = sum(diag(A)) # trace of hat matrix
rss = crossprod(y - f_hat) # residual sum of squares
gcv = n*rss / (n - tr_A)^2 # compute GCV score
return(list(f_hat = f_hat, omega_tilde = omega_tilde, gcv = gcv, rss=rss))
}Using the ALS dataset, we will implement a direct search for the GCV optimal ridge regularization parameter \(\lambda\) to predict the clinical outcome ALSFRS_slope. We will use the R package furrr and request multiple cores (\(workers = 8\)) in the call to furrr::future_map(). Using multicore processing will significantly expedite this massive search calculation.
We use LASSO (alpha = 1) regularization on the same ALS test case with
a Laplacian reproducing kernel operator, rk()
\[\kappa(s,t) = rk(s,t)= e^{-\left (\frac{\lVert s - t\rVert}{\sigma}\right )}\ .\]
rk_Laplace <- function(s, t, sigma=1) {
if (sigma <= 0) sigma=1 # avoid singularities
init_len = length(s)
rk_Lap = 0
for (i in 1:init_len) { rk_Lap = (s[i] - t[i])^2 + rk_Lap }
rk_Lap = exp(-(sqrt(rk_Lap))/sigma)
return(rk_Lap)
} Next we fit the Ridge model and plot the trajectories of the parameter estimates across
\(\lambda\).
lambda = 10^seq(-6, 2, by = 0.5)
library(furrr)
# For classical sequential parameter search use:
# gcv_search = map(lambda, function(lam) ridge(X, y, lam))
# For multicore execution use:
plan(multisession, workers = 8)
gcv_search <- future_map(lambda, function(lam) ridge(X, y, lam,rk_Laplace))
V = map_dbl(gcv_search, function(x) x$gcv)
ridge_coefs = map_df(gcv_search, function(x) {
data.frame(value = x$omega_tilde, coef =as.character(1:2223))
},
.id = 'iter') %>%
mutate(lambda = lambda[as.integer(iter)])library(patchwork)
library(plotly)
# gcv_plot = qplot(lambda, V, geom='line', main='GCV score (log scale)', ylab='GCV') +
# scale_x_log10()
plot_ly(x=lambda, y=V, type="scatter", mode="lines") %>%
layout(title="(Ridge) GCV vs. Log10(Lambda)",
xaxis = list(title="log10(Lambda)", type = "log"))# beta_plot = ridge_coefs %>%
# ggplot(aes(x = lambda, y = value, color = coef)) +
# geom_line() +
# scale_x_log10() +
# scico::scale_color_scico_d(end = 0.8) +
# labs(title = 'Model Effects (Betas) Across Lambda')
# ggplotly(gcv_plot + beta_plot)
plot_ly(data=ridge_coefs[as.numeric(ridge_coefs$coef)<20,], x=~lambda, y=~value,
type="scatter", mode="lines", color=~coef, name=~coef) %>%
layout(title="(Ridge Estimates) Effects (Beta's) vs. log10(Lambda)",
xaxis = list(title="", type = "log"), # title="log10(Lambda)",
legend = list(title="Predicting Covariates", orientation = 'h'))The plot allows manual zoom in to confirm that instances (case observations) \(4\) and \(18\) appear to be most dominating in estimating the case-contribution weights in the linear regularized model, among the first \(20\) participants. By zooming in or deselecting these two instances \(\{4,18\}\), we can explore the canonical regularization path trajectories.
The following example uses small numbers to reduce the extreme computational burden and expedite the computational estimates in compiling the HTML output from the Rmd source. In practice, larger numbers will provide more precise estimates. The resulting RSS are reported in tables.
# Matern kernel function
maternkernel <- function(sigma,nu,rho,x){
scaled_x = sqrt(2*nu)*x/rho
return (sapply(scaled_x,besselk_nu(nu))*scaled_x**nu*sigma**2*2**(1-nu)/gamma(nu))
}
rk_Matern <- function(s, t, nu=0.5) {
# if (sigma <= 0) sigma=1 # avoid singularities
init_len = length(s)
rk_dist = 0
for (i in 1:init_len) { rk_dist = (s[i] - t[i])^2 + rk_dist }
rk_Mat = maternkernel(1,nu,1, rk_dist+1e-6)
return(rk_Mat)
}
maternkernel_nu_mod <- function(nu) {
return(function(s, t) {
# if (sigma <= 0) sigma=1 # avoid singularities
init_len = length(s)
rk_dist = 0
for (i in 1:init_len) { rk_dist = (s[i] - t[i])^2 + rk_dist }
rk_Mat = maternkernel(1,nu,1, sqrt(rk_dist)+1e-6)
return(rk_Mat)
} )
}rk_Matern <- function(s, t, nu=0.5) {
# if (sigma <= 0) sigma=1 # avoid singularities
init_len = length(s)
rk_dist = 0
for (i in 1:init_len) { rk_dist = (s[i] - t[i])^2 + rk_dist }
rk_Mat = maternkernel(1,nu,1, sqrt(rk_dist)+1e-6)
return(rk_Mat)
} poly_d_mod <- function(d) {
return(function(s, t) {
init_len = length(s)
rk = 1
for (i in 1:init_len) { rk = s[i]*t[i] + rk }
return(rk**d)
} )
}Next, we will expand the previous standard Ridge model by using a RKHS and specifying alternative kernels in the estimation of the model coefficients.
For \(\lambda\in \{0.01, 0.1, 1, 10\}\).
#Hyperparameter set 1 - Linear Kernel
library(foreach)
library(doParallel)
registerDoParallel(8)
# fitted_ridge1<- ridge(X,y,lambda=10)
# fitted_ridge2<- ridge(X,y,lambda=1)
# fitted_ridge3<- ridge(X,y,lambda=0.1)
# fitted_ridge4<- ridge(X,y,lambda=0.01)
lambda <- c(0.001, 0.01, 0.05, 0.1, 0.5, 1, 2, 10)
fitted_ridge <- vector()
fitted_ridge <- foreach (param = 1:length(lambda), .combine = 'c') %dopar% {
ridge(X,y,lambda=lambda[param])
}
stopImplicitCluster()df <- as.data.frame(fitted_ridge)
df_tbl <- rbind(c(df$gcv[1], df$gcv.1[1], df$gcv.2[1], df$gcv.3[1],
df$gcv.4[1], df$gcv.5[1], df$gcv.6[1], df$gcv.7[1]),
c(df$rss[1], df$rss.1[1], df$rss.2[1], df$rss.3[1]))
colnames(df_tbl) <- c(paste0("lambda=", lambda[1]), paste0("lambda=", lambda[2]),
paste0("lambda=", lambda[3]), paste0("lambda=", lambda[4]),
paste0("lambda=", lambda[5]), paste0("lambda=", lambda[6]),
paste0("lambda=", lambda[7]), paste0("lambda=", lambda[8]))
rownames(df_tbl) <- c("GCV", "RSS")
library(DT)
datatable(df_tbl, caption="Linear Ridge Kernel (parameter estimates)")
# fitted_ridge1$rss
# fitted_ridge2$rss
# fitted_ridge3$rss
# fitted_ridge4$rssFor \(\lambda=0.1\) and \(\nu\in \{0.01, 0.1, 0.5, 1, 1.5, 2, 2.5, 10\}\).
# Hyperparameter set 2 - Matern kernel
# matern1 <- ridge(X,y,lambda=0.1,kernel=rk_Matern)
# matern2 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(1))
# matern3 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(0.01))
# matern4 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(1.5))
# matern5 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(2.5))
# matern6 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(10))
# matern7 <- ridge(X,y,lambda=0.1,kernel=maternkernel_nu_mod(0.5))
registerDoParallel(8)
nu <- c(0.01, 0.1, 0.5, 1, 1.5, 2, 2.5, 10)
fitted_ridgeMatern <- vector()
fitted_ridgeMatern <-
foreach (param = 1:(length(nu)+1), .combine = 'c') %dopar% {
besselk_nu <- function(nu) {
return(function(x) { return(GPBayes::BesselK(nu, x)) }) }
if (param==1) ridge(X, y, lambda=0.1, kernel=rk_Matern)
else ridge(X, y, lambda=0.1, kernel=maternkernel_nu_mod(nu[param-1]))
}
stopImplicitCluster()
df <- as.data.frame(fitted_ridgeMatern)
df_tbl <- rbind(c(df$gcv[1], df$gcv.1[1], df$gcv.2[1], df$gcv.3[1], df$gcv.4[1],
df$gcv.5[1], df$gcv.6[1], df$gcv.7[1], df$gcv.8[1], df$gcv.9[1]),
c(df$rss[1], df$rss.1[1], df$rss.2[1], df$rss.3[1], df$rss.4[1],
df$rss.5[1], df$rss.6[1], df$rss.7[1], df$rss.8[1], df$rss.9[1]))
colnames(df_tbl) <- c(paste0("nu=", 0.5), paste0("nu=", nu[1]), paste0("nu=", nu[2]),
paste0("nu=", nu[3]), paste0("nu=", nu[4]),
paste0("nu=", nu[5]), paste0("nu=", nu[6]),
paste0("nu=", nu[7]), paste0("nu=", nu[8]))
rownames(df_tbl) <- c("GCV", "RSS")
library(DT)
datatable(df_tbl, caption="Matern Ridge Kernel (parameter estimates)")We can also use any of the other kernels (e.g., polynomial) to fit RKHS Ridge models and compare their performance across different kernels, and hyperparameters (e.g., \(\nu\), \(d\)).
#Hyperparameter set 3 - Polynomial Kernel kernel
# poly1 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(3/2))
# poly2 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(1))
# poly3 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(0.9))
# poly4 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(0.7))
# poly5 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(0.5))
# poly6 <- ridge(X,y,lambda=0.1,kernel=poly_d_mod(0.3))
registerDoParallel(7)
poly_d <- c(0.3, 0.5, 0.7, 0.9, 1.5, 2.5, 10)
fitted_ridgePoly <- vector()
fitted_ridgePoly <-
foreach (param = 1:length(poly_d), .combine = 'c') %dopar% {
ridge(X, y, lambda=0.1, kernel=poly_d_mod(poly_d[param]))
}
stopImplicitCluster()
df <- as.data.frame(fitted_ridgePoly)
df_tbl <- rbind(c(df$gcv[1], df$gcv.1[1], df$gcv.2[1], df$gcv.3[1], df$gcv.4[1],
df$gcv.5[1], df$gcv.6[1]),
c(df$rss[1], df$rss.1[1], df$rss.2[1], df$rss.3[1], df$rss.4[1],
df$rss.5[1], df$rss.6[1]))
colnames(df_tbl) <- c(paste0("d=", poly_d[1]), paste0("d=", poly_d[2]),
paste0("d=", poly_d[3]), paste0("d=", poly_d[4]),
paste0("d=", poly_d[5]), paste0("d=", poly_d[6]),
paste0("d=", poly_d[7]))
rownames(df_tbl) <- c("GCV", "RSS")
library(DT)
datatable(df_tbl, caption="Polynomial Ridge Kernel (parameter estimates)")Finally, we can select the best model corresponding to the optimal regularization parameter \(\lambda=\lambda_{[which.min(V)]}\) associated with minimization of the GCV search and obtain the corresponding effect size estimates \(\beta\)’s, for example, by refitting the Ridge model with a linear kernel.
plan(multisession, workers = 8)
X_augmented <- cbind(1, X) # augmented design matrix
gcv_search <- future_map(lambda, function(lam) ridge(X_augmented, y, lam))
V = map_dbl(gcv_search, function(x) x$gcv)
ridge_coefs = map_df(gcv_search, function(x) {
data.frame(value = x$omega_tilde, coef =as.character(1:2223))
},
.id = 'iter') %>% mutate(lambda = lambda[as.integer(iter)])
fit_ridge <- ridge(X_augmented, y, lambda[which.min(V)]) # fit optimal model
omega_tilde <- fit_ridge$omega_tilde
beta_hat <- crossprod(omega_tilde, X) # slope and noise term coefficientsNext, we will compare our manual ridge-regularized linear model parameter estimation strategy (using RKHS) against the classical glmnet() gold standard. Recall that setting the glmnet hyperparameter alpha = 0 and alpha = 1 correspond to ridge (\(L_2\)) and LASSO (\(L_1\)) regression, respectively.
library(glmnet)
# c(beta_0, beta_hat)
fit_glmnet = glmnet(X, y, alpha = 0, lambda=lambda, standardize = FALSE)The following table contrasts the manual RKHS-Ridge and the automated
GLM-Net Ridge regression effect size estimates for all \(5\) covariate
predictors of the outcome ALSFRS_slope.
Theoretically, the manual estimates are not expected to match
the GLM-Net Ridge estimates, as \(X^TX\neq I\).
kable_df = function(..., digits =3, full_width = FALSE)
kable(..., digits=digits) %>% kable_styling(full_width = full_width)
rbind(RKHS_Ridge = c(beta_hat), # -1 to skip over the intercept estimation
GLMNet_Ridge = (coef(fit_glmnet)[, which.max(fit_glmnet$dev.ratio)])[-1]) %>%
kable_df()| Age_mean | Albumin_max | Albumin_median | Albumin_min | Albumin_range | |
|---|---|---|---|---|---|
| RKHS_Ridge | 0.042 | 1.277 | -0.313 | -0.246 | -4.674 |
| GLMNet_Ridge | 0.042 | 1.304 | -0.319 | -0.261 | -4.747 |
Again, we will define a manual lasso() method de novo using RKHS modeling,
the Gram matrix, and fitting a LASSO regularized linear model for estimating the parameters (effect sizes). To see the theoretical derivation for kernelized LASSO regression (See appendix 5.4).
There is a distinction of shrinking on observations and shrinking on features. We run some experiments showing a possible way of implementing kernelized LASSO with shrinkage on the features or on the instances. Then, we will compare the results with the classical glmnet approach for efficient elastic net regularized linear modeling.
To demonstrate participant selection, i.e., instance or case shrinkage, we’ll use the following linear kernel implementation.
library(MASS)
# observation shrinkage
lasso_observation <- function(X, y, lambda=0.1, kernel=rk) {
K = gram(X,kernel)
# assumed X to be n*d where observations stacked on the rows,
# GramMatrix (nxn) XTX (already phi^Tphi)
n = length(y)
#Q = cbind(1, GramMatrix) no need for intercept as the kernel already captures this
#S = rbind(0, cbind(0, GramMatrix))
alpha_hat = inverse(K) %*% y
for (j in 1:n) { # The explicit formula for LASSO loss
alpha_hat[j] = alpha_hat[j]*max(0,1-length( alpha_hat)*lambda/abs( alpha_hat[j]))
}
M = K%*%K + lambda*ginv(diag(as.list(abs(alpha_hat)))) # X^TX+lambda W-
M_inv = inverse(M) # invert M
f_hat = K %*% alpha_hat
A = K %*% M_inv %*% K # Hat matrix X(X^XX+lambda W-)^_1 X^T
tr_A = sum(diag(A)) # trace of hat matrix
rss = crossprod(y - f_hat) # residual sum of squares
gcv = n*rss / (n - tr_A)^2 # compute GCV score
return(list(f_hat = f_hat, alpha_hat = alpha_hat, gcv = gcv, rss=rss))
}The prediction procedure is essentially \[\hat{y}=K S_{\lambda}(K^{-1}y)\ .\] Notice that when \(\lambda=0\), this degenerates into a perfect reconstruction. As a benchmark, the resulting RSS and GCV values are \(34.57\) and \(11.71\), respectively.
# when lambda is 0
fitted_lambda0 <- lasso_observation(X, y, 0,rk_Laplace)
fitted_lambda0$gcv
fitted_lambda0$rssNext, we run over the grid space of the \(\lambda\) parameter.
lambda = 10^seq(-9, 2, by = 0.5)
library(furrr)
# For classical sequential parameter search use:
# gcv_search = map(lambda, function(lam) ridge(X, y, lam))
# For multicore execution use:
plan(multisession, workers = 12)
gcv_search <- future_map(lambda, function(lam) lasso_observation(X, y, lam,rk_Laplace))
V_lasso = map_dbl(gcv_search, function(x) x$gcv)
lasso_coefs = map_df(gcv_search, function(x) {
data.frame(value = x$alpha_hat, coef =as.character(1:2223))
},
.id = 'iter') %>% mutate(lambda = lambda[as.integer(iter)])We observe that initially, the GCV score decreases, since \(GCV(lambda=0)= 11.71\)), and then alternates between increases and decreases.
library(patchwork)
library(plotly)
RSS_lasso = map_dbl(gcv_search, function(x) x$rss)
plot_ly(x=lambda, y=RSS_lasso, type="scatter", mode="lines") %>%
layout(title="(LASSO) RSS (Observation) vs. Log10(Lambda)",
xaxis = list(title="log10(Lambda)", type = "log"))# gcv_plot = qplot(lambda, V, geom='line', main='GCV score (log scale)', ylab='GCV') +
# scale_x_log10()
plot_ly(x=lambda, y=V_lasso, type="scatter", mode="lines") %>%
layout(title="(LASSO) GCV (Observation) vs. Log10(Lambda)",
xaxis = list(title="log10(Lambda)", type = "log"))# beta_plot = ridge_coefs %>%
# ggplot(aes(x = lambda, y = value, color = coef)) +
# geom_line() +
# scale_x_log10() +
# scico::scale_color_scico_d(end = 0.8) +
# labs(title = 'Model Effects (Betas) Across Lambda')
# ggplotly(gcv_plot + beta_plot)
plot_ly(data=lasso_coefs[as.numeric(lasso_coefs$coef)<20,], x=~lambda, y=~value,
type="scatter", mode="lines", color=~coef, name=~coef) %>%
layout(title="LASSO Effect Estimates (Beta's) vs. log10(Lambda)",
xaxis = list(title="", type = "log"), # title="log10(Lambda)",
legend = list(title="Predicting Covariates", orientation = 'h'))index <- 1
meta_matrix <- matrix(0,nrow=2223,ncol=length(lambda)) # 2223 is the number of rows
for (i in lambda){
meta_matrix[,index] <- (lasso_coefs[lasso_coefs$lambda==i,]$value ==0)
index <- index + 1
}
library(DT)
# image(t(meta_matrix),xlab="Lambda values", yaxt="n",xaxt="n")
# axis(3, at=seq(0,1, length=length(lambda)), labels=as.character(signif(lambda,2)))
plot_ly(z=t(meta_matrix), type="heatmap") %>%
layout(title="Lambda values (x-axis) vs. LASSO Coefficient Estimates (y-axis)",
xaxis = list(tickvals = seq(0, dim(meta_matrix)[1], length=length(lambda)),
ticktext = c(as.character(signif(lambda,2)))))The large lambda values (on the right) indicate trivial effect size estimates, whereas the smaller lambda values (on the left) signify nonzero effect size estimates for the importance of the observation instance’s (cases) as contributors to the LASSO model prediction. Of course, due to the shrinkage nature for LASSO, as lambda increases, all observation effect sizes are discounted and thresholded to \(0\).
Selection of salient features (variable selection) is independent from the prior search for influential cases (observations) that impact the LASSO model prediction.
Now we focus on the traditional feature selection using LASSO.
kernelized_y <- function(X,y, rkfunc = rk) { # compute the `crossprod` using the specified RK
apply(X, 1, function(Row)
rkfunc(Row, y) # specifies the Reproducing Kernel
)
}library(MASS)
# observation shrinkage
lasso_feature <- function(X, y, lambda=0.1, kernel=rk) {
K = gram(t(X),kernel) #assumed X to be n*d where observations stacked on the rows, GramMatrix (nxn) XTX (already phi^Tphi)
n = ncol(X)
#Q = cbind(1, GramMatrix) no need for intercept as the kernel already captures this
#S = rbind(0, cbind(0, GramMatrix))
kernel_yval <-kernelized_y(t(X),y,kernel)
beta_hat = inverse(K) %*% kernel_yval
for (j in 1:n) { # The explicit formula for LASSO loss
beta_hat[j] = beta_hat[j]*max(0,1-length(beta_hat )*lambda/abs(beta_hat[j]))
}
M = K%*%K + lambda*ginv(diag(as.list(abs(beta_hat)))) # X^TX+lambda W-
M_inv = inverse(M) # invert M
f_hat = K %*% beta_hat
A = K %*% M_inv %*% K # Hat matrix X(X^XX+lambda W-)^_1 X^T
tr_A = sum(diag(A)) # trace of hat matrix
rss = crossprod(kernel_yval - f_hat) # residual sum of squares
gcv = n*rss / (n - tr_A)^2 # compute GCV score
return(list(f_hat = f_hat, beta_hat = beta_hat, gcv = gcv, rss=rss))
}Let’s examine the GCV and RSS performance metrics.
Fit the LASSO feature-selection model over the \(\lambda\) space partition (grid).
lambda = 10^seq(-9, 2, by = 0.5)
library(furrr)
# For classical sequential parameter search use:
# gcv_search = map(lambda, function(lam) ridge(X, y, lam))
# For multicore execution use:
plan(multisession, workers = 12)
gcv_search <- future_map(lambda, function(lam) lasso_feature(X_augmented, y, lam))
V_lasso = map_dbl(gcv_search, function(x) x$gcv)
lasso_coefs = map_df(gcv_search, function(x) {
data.frame(value = x$beta_hat, coef =as.character(1:6))
},
.id = 'iter') %>% mutate(lambda = lambda[as.integer(iter)])Plot the resulting effect-size (beta) estimate trajectories.
library(patchwork)
library(plotly)
RSS_lasso = map_dbl(gcv_search, function(x) x$rss)
plot_ly(x=lambda, y=RSS_lasso, type="scatter", mode="lines") %>%
layout(title="(LASSO) RSS (Feature) vs. Log10(Lambda)",
xaxis = list(title="log10(Lambda)", type = "log"))# gcv_plot = qplot(lambda, V, geom='line', main='GCV score (log scale)', ylab='GCV') +
# scale_x_log10()
plot_ly(x=lambda, y=V_lasso, type="scatter", mode="lines") %>%
layout(title="(LASSO) GCV (Feature) vs. Log10(Lambda)",
xaxis = list(title="log10(Lambda)", type = "log"))# beta_plot = ridge_coefs %>%
# ggplot(aes(x = lambda, y = value, color = coef)) +
# geom_line() +
# scale_x_log10() +
# scico::scale_color_scico_d(end = 0.8) +
# labs(title = 'Model Effects (Betas) Across Lambda')
# ggplotly(gcv_plot + beta_plot)
plot_ly(data=lasso_coefs, x=~lambda, y=~value,
type="scatter", mode="lines", color=~coef, name=~coef) %>%
layout(title="(LASSO Estimates) Effects (Beta's) vs. log10(Lambda)",
xaxis = list(title="", type = "log"), # title="log10(Lambda)",
legend = list(title="Predicting Covariates", orientation = 'h'))index <- 1
meta_matrix <- matrix(0,nrow=6,ncol=length(lambda)) # 2223 is the number of columns
for (i in lambda){
meta_matrix[,index] <- (lasso_coefs[lasso_coefs$lambda==i,]$value ==0)
index <- index + 1
}
library(DT)
# image(t(meta_matrix),xlab="Lambda values", yaxt="n",xaxt="n")
# axis(3, at=seq(0,1, length=length(lambda)), labels=as.character(signif(lambda,2)))
plot_ly(z=t(meta_matrix), type="heatmap") %>%
layout(title="Lambda values (x-axis) vs. LASSO Coefficient Estimates (y-axis)",
xaxis = list(tickvals = seq(0, dim(meta_matrix)[1], length=length(lambda)),
ticktext = c(as.character(signif(lambda,2)))))Direct comparisons between different RSS measures across different kernels is challenging, as the corresponding feature maps \(\phi(y)\) are different and hence, the RKHS embedding is in different space scales. For instance, the scale of the Laplace kernel feature map may be very small compared to the larger scale of the feature map associated with the linear kernel. There is no uniform criterion to interpreting performance metrics such as RSS across different RKHS embedding spaces.
As a sanity check at the end, we can select the best model corresponding to the
optimal regularization parameter \(\lambda=\lambda_{[which.min(V)]}\) associated
with minimization of the GCV search and obtain the corresponding effect size
estimates \(\beta\)’s. This will allow us to contrast our implementation of
LASSO kernel feature selection shrinkage against the linear dot product kernel,
which roughly matches the glmnet() approach. Again, these will not be
identical as \(X^TX\neq I\).
Let’s compare the manual LASSO-regularized linear model parameter estimation strategy (using RKHS) against the classical glmnet(), as a gold standard. Recall that setting the glmnet hyperparameter alpha = 1 corresponds to LASSO (\(L_1\)) regression.
library(glmnet)
# c(beta_0, beta_hat)
fit_glmnet = glmnet(X, y, alpha = 1, lambda=lambda, standardize = FALSE)The following table compares the manual RKHS-LASSO and the automated
GLM-Net LASSO regression effect size estimates for all \(5\) covariate
predictors of the outcome ALSFRS_slope.
kable_df = function(..., digits =3, full_width = FALSE)
kable(..., digits=digits) %>% kable_styling(full_width = full_width)
rbind(RKHS_LASSO = c(beta_hat),
GLMNet_LASSO = coef(fit_glmnet)[, which.max(fit_glmnet$dev.ratio)]) %>%
kable_df()| (Intercept) | Age_mean | Albumin_max | Albumin_median | Albumin_min | Albumin_range | |
|---|---|---|---|---|---|---|
| RKHS_LASSO | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | -2.861 |
| GLMNet_LASSO | -0.518 | 0.042 | 1.308 | -0.319 | -0.264 | -4.758 |
Compare these LASSO (effect-size) \(\beta\) estimates against their Ridge counterparts, which we computed earlier for the same ALS data.
For this example we’ll simulate a spirals dataset representing 2D point cloud data of two spirals corrupted with Gaussian noise using the functions kernlab::spirals()
and mlbench::mlbench.spirals().
### Binary Cross-Entropy (BCE) loss function
The binary cross-entropy (BCE) loss function, also called Log Loss, or logistic loss can be used for binary outcome predictions
\[BCE\equiv L_{\log}(y,p) = -\left ( y\ \log(p) + (1-y)\ \log(1-p) \right )\ , \]
where \(y\) is the true binary label (e.g., \(0\) or \(1\)) and \(p\) is the model-based probability estimate that \(0\leq p=Prob(y=1)\leq 1\). A more general Balanced Log Loss includes weight coefficients \(w(y_i)\)
\[BalancedBCE\equiv = -\frac{1}{n}\sum_{i=1}^n { w(y_i)\left ( y_i\ \log(p(y_i)) + (1-y(p_i))\ \log(1-p(y_i)) \right )}\ , \]
where \(n\) is the total number of samples, \(y_i\) is the true binary label of the \(i^{th}\) sample (\(0\) or \(1\)), \(p(y_i)\) is the predicted probability of the \(i^{th}\) sample being of class \(1\), and the weights \(w(y_i)\) for class \(y_i\) are appropriately chosen, e.g., the inverse of the class frequencies.
To compute the Gram matrix, we need to first define the spline reproducing kernel function rk_spline() implementing the special case of a piecewise cubic polynomial
\[\kappa(x, y) = \prod_{i=1}^n\left (1+ x_i y_i (1+\min(x_i,y_i))- \frac{(x_i+y_i)\min^2(x_i,y_i)}{2} +\frac{\min^3(x_i,y_i)}{3}\right )\ .\]
Next we define the smoothing_spline() function using the Gram matrix and the
RKHS rk_spline() function.
# May need to use BCE loss function (instead of RSS)
# binary_cross_entropy(pred, data.valid,loss.unit=c("individuals","L2 units"), y, L2.unit)
# https://rdrr.io/cran/autoMrP/src/R/utils.R
###########################################################################
# binary cross-entropy ----------------------------------------------------
###########################################################################
binary_cross_entropy <- function(dataProb, trueLabels, positiveClassLabel="1",
w_positive=0.5, w_negative=0.5) {
# positiveClassLabel is the positive class
# Extract true labels and predicted probabilities
truth <- as.numeric(trueLabels == positiveClassLabel)
prob <- dataProb # Probabilities for positive class
# Calculate the BCE balanced log loss
BCE <- -mean(
(truth * log(prob) * w_positive) + ((1 - truth) * log(1 - prob) * w_negative)
)
return(BCE)
}
# Normalize all probability estimates to [0-e, 1-e]
probNormalize <- function(p) {
epsil <- 1e-4 # avoid singularities
m <- min(p)
M <- max(p)
return ( (p-m+epsil)/(M-m+2*epsil) )
}
# # ################# BCE test
# # # smoothing_spline(x, y, lam)
# X = x
# y = as.numeric(labels)
#
# GramMatrix = gram(X, rkfunc = rk_spline) # GramMatrix matrix (nxn)
#
# n = length(y)
# J = cbind(1, X) # matrix with a basis for the null space of the penalty
# Q = cbind(J, GramMatrix) # design matrix
# m = ncol(J) # dimension of the null space of the penalty
#
# S = matrix(0, n + m, n + m) # initialize S
# S[(m + 1):(n + m), (m + 1):(n + m)] = GramMatrix # non-zero part of S
#
# M = crossprod(Q) + lambda[1]*S
# M_inv = inverse(M) # inverse of M
#
# gamma_hat = crossprod(M_inv, crossprod(Q, y))
# f_hat = Q %*% gamma_hat
#
# A = Q %*% M_inv %*% t(Q)
# tr_A = sum(diag(A)) # trace of hat matrix
#
# # rss = crossprod(y - f_hat) # BCE loss
# # pred <- as.factor(ifelse (f_hat > 1.5, "1", "2"))
# pred <- probNormalize(f_hat - 1)
# BCE <- binary_cross_entropy(dataProb=pred[, 1], trueLabels=as.character(labels),
# positiveClassLabel="1")
#########################################
# Numerical output (regression spline)
#########################################
smoothing_spline <- function(X, y, lambda) {
GramMatrix = gram(X, rkfunc = rk_spline) # GramMatrix matrix (nxn)
n = length(y)
J = cbind(1, X) # matrix with a basis for the null space of the penalty
Q = cbind(J, GramMatrix) # design matrix
m = ncol(J) # dimension of the null space of the penalty
S = matrix(0, n + m, n + m) # initialize S
S[(m + 1):(n + m), (m + 1):(n + m)] = GramMatrix # non-zero part of S
M = crossprod(Q) + lambda*S
M_inv = inverse(M) # inverse of M
gamma_hat = crossprod(M_inv, crossprod(Q, y))
f_hat = Q %*% gamma_hat
A = Q %*% M_inv %*% t(Q)
tr_A = sum(diag(A)) # trace of hat matrix
rss = crossprod(y - f_hat) # residual sum of squares (RSS)
gcv = n * rss/(n - tr_A)^2 # compute GCV score
return(list(f_hat = f_hat, gamma_hat = gamma_hat, gcv = gcv))
}
#########################################
# Categorical output (binary cross-entropy spline)
#########################################
smoothing_BCE_spline <- function(X, y, lambda) {
GramMatrix = gram(X, rkfunc = rk_spline) # GramMatrix matrix (nxn)
n = length(y)
J = cbind(1, X) # matrix with a basis for the null space of the penalty
Q = cbind(J, GramMatrix) # design matrix
m = ncol(J) # dimension of the null space of the penalty
S = matrix(0, n + m, n + m) # initialize S
S[(m + 1):(n + m), (m + 1):(n + m)] = GramMatrix # non-zero part of S
M = crossprod(Q) + lambda*S
M_inv = inverse(M) # inverse of M
gamma_hat = crossprod(M_inv, crossprod(Q, y))
f_hat = Q %*% gamma_hat
A = Q %*% M_inv %*% t(Q)
tr_A = sum(diag(A)) # trace of hat matrix
# rss = crossprod(y - f_hat) # BCE loss
# pred <- as.factor(ifelse (f_hat > 1.5, "1", "2"))
# convert fitted values to probabilities in [0, 1] for BCE loss estimation
pred <- probNormalize(f_hat - 1) # 1 <= f_hat <=2
BCE <- binary_cross_entropy(dataProb=pred[, 1], trueLabels=as.character(labels),
positiveClassLabel="1")
return(list(f_hat = f_hat, gamma_hat = gamma_hat, BCE = BCE))
}Now we can fit the RKHS spline model for multiple \(\lambda\) values and optimize for the GCV metric to select the optimal \(\lambda\) parameter.
lambda = 10^seq(-6, 0, by = 0.5)
# Slow # gcv_search = map(lambda, function(lam) smoothing_spline(x, y, lam))
# For quicker multi core execution
plan(multisession, workers = 8)
gcv_search = map(lambda, function(lam) smoothing_spline(x, y, lam))
V = map_dbl(gcv_search, function(x) x$gcv)
## BCE Metric
plan(multisession, workers = 8)
BCE_search = map(lambda, function(lam) smoothing_BCE_spline(x, y, lam))
V_BCE = map_dbl(BCE_search, function(x) x$BCE)Plot of GCV across search space of \(\lambda\) values.
Again, we can compare and contrast the manual spline RKHS model fit estimation
against a standard model, in this case, a generalized additive model (GAM)
estimated using function mgcv::gam(). See the gam model specification format.
The graph below shows the raw spirals simulated data (as scatter points color coded by their true arc-type class labels), and three different models - Spline-RKHS, GAM, and Binary classifier. The shapes and colors of the data-points and the predictions of the \(3\) models (also rendered as scatter to avoid occlusions/overlaps) represent the true class labels (shapes) and the model-fit (colors). Note how well the binary-GAM model discriminates between the two spiral arches, correctly identifying the majority of the yellow-arch points (high positive values, red-triangle symbols) and black-arch points (low negative values, red-circle symbols).
library(mgcv)
# https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/gam.models.html
fit_rk = smoothing_spline(x, y, lambda[which.min(V)]) # fit optimal model
fit_rk_bin = smoothing_BCE_spline(x, y, lambda[which.min(V_BCE)]) # fit optimal model
gamSpline = gam(y ~ s(x,bs="cr",k=100)) # s(x) = smooth function of x
print(paste0("y ~ s(x) GAM Spline model Log-Likelihood = ", round(logLik(gamSpline), 4)))## [1] "y ~ s(x) GAM Spline model Log-Likelihood = -339.812"gamBin <- gam(labels ~ s(x,y), family = 'binomial')
print(paste0("labels ~ s(x,y) GAM Binary model Log-Likelihood = ", round(logLik(gamBin), 4)))## [1] "labels ~ s(x,y) GAM Binary model Log-Likelihood = -135.8485"df <- as.data.frame(cbind(X=spirals$x[ , 1], Y=spirals$x[ , 2],
class=as.factor(spirals$classes)))
# fit_plot = df %>%
# mutate(fit_rk = fit_rk$f_hat[, 1],
# fit_gam = fitted(fit_gam)) %>%
# arrange(X) %>%
# pivot_longer(-c(X, Y), names_to = 'fit', values_to = 'value') %>%
# ggplot(aes(X, Y)) +
# geom_point(color = '#FF55001A') +
# geom_line(aes(y = value, color = fit)) +
# scico::scale_color_scico_d(palette = 'hawaii', begin = .2, end = .8)
# ggplotly(gcv_plot + fit_plot)
# identify the 2 true spirals
# df_reformat$class <- rep(spirals$classes, 2)
# table(fit_rk_bin$f_hat[, 1] > 1.499)
# FALSE TRUE
# 152 148
df_reformat <- df %>%
mutate(RKHS_mod=fit_rk$f_hat[, 1],
RKHS_Spline_Binary=2*as.numeric(fit_rk_bin$f_hat[, 1] > 1.499)-1,
GAM_Spline=fitted(gamSpline), GAM_Binary=3*fitted(gamBin)-1.5) %>%
arrange(X) %>%
pivot_longer(-c(X, Y, class), names_to = 'fit', values_to = 'value')
df_reformat %>%
plot_ly(x=~X, y=~Y, type="scatter", mode="markers", color=~class,
# symbol=~class, symbols=as.character(seq_along(unique(df_reformat$class))),
name=~paste0("True-Class=", class)) %>% hide_colorbar() %>%
add_trace(x=~X, y=~value, color=~fit, colors = c("red", "blue", "green", "purple"),
symbol=~class, symbols=as.character(seq_along(unique(df_reformat$class))),
name=~paste0(" ", fit), mode="markers",
hovertemplate = paste('(Model)')) %>%
# colors = RColorBrewer::brewer.pal(3, "Set2")[1:2]) %>%
layout(title="(Binary & Regressive) GAM vs. RKHS Models of Spirals Data\n
Color=Model Fit, Shape=True-Class",
xaxis = list(title="X", showline=TRUE), yaxis = list(title="Y"),
legend = list(title="Objects", orientation = 'h'), hovermode = 'x')In this section we will connect spacekime analytics with Reproducing Kernel Hilbert Spaces (RKHS). The basic idea is that repeated measurements of longitudinal processes can be viewed as time-varying probability distributions. The explicit connection is that the kime-phase distribution \(\Phi_{[-\pi,\pi)}^{(t)}\) represents the distribution of the repeated measurements at a fixed spatiotemporal location \(t\). As time (kime-magnitude) and space evolve, these phase distributions change and their dynamics can be modeled using RKHS. Recall the following example shown in TCIU Chapter 3 appendix (Radon-Nikodym Derivatives, Kime Measures, and Kime Operator.
In this simulation, the \(3\) spatial dimensions \((x,y,z)\in\mathbb{R}^3\) are compressed into \(1D\) along the vertical axis \(z\in\mathbb{R}^1\), the radial displacement represents the time dynamics, \(t\in \mathbb{R}^+\), and the angular scatters of three different processes, representing repeated measurements from \(3\) different circular distributions at a fixed spatiotemporal location, are shown in different colors. Mean-dynamics across time of the three different time-varying distributions are shown as smooth curves color-coded to reflect the color of their corresponding phase distributions, samples, and sampling-distributions.
library(animation)
library(circular)
library(plotly)
epsilon <- 0.1
sampleSize <- 1000 # total number of phases to sample for 3 different processes (x, y, z)
sizePerTime <- 100 # number of phases to use for each fixed time (must divide sampleSize)
circleUniformPhi <- seq(from=-pi, to=pi, length.out=sizePerTime)
oopt = ani.options(interval = 0.2)
set.seed(1234)
# sample the the kime-phases for all 3 different processes and the r time points
x <- rvonmises(n=sampleSize, mu=circular(pi/5), kappa=3)
y <- rvonmises(n=sampleSize, mu=circular(-pi/3), kappa=5)
z <- rvonmises(n=sampleSize, mu=circular(0), kappa=10)
r <- seq(from=1, to=sampleSize/sizePerTime, length.out=10)
# Define a function that renormalizes the kime-phase to [-pi, pi)
pheRenormalize <- function (x) {
out <- ifelse(as.numeric(x) <= pi, as.numeric(x)+pi, as.numeric(x)-pi)
return (out)
}
# transform Von Mises samples from [0, 2*pi) to [-pi, pi)
x <- pheRenormalize(x)
y <- pheRenormalize(y)
z <- pheRenormalize(z)
# vectorize the samples
vectorX = as.vector(x)
vectorY = as.vector(y)
vectorZ = as.vector(z)
# Starting phases, set the first phase index=1
plotX = c(vectorX[1])
plotY = c(vectorY[1])
plotZ = c(vectorZ[1])
# pl_list <- list()
pl_scene <- plot_ly(type='scatter3d', mode="markers")
plotX <- list()
plotY <- list()
plotZ <- list()
plotX_df <- list() # need separate data frames to store all time foliations
plotY_df <- list()
plotZ_df <- list()
for (t in 1:length(r)) { # loop over time
# loop over kime-phases
plotX[[t]] <- as.numeric(x[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
plotY[[t]] <- as.numeric(y[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
plotZ[[t]] <- as.numeric(z[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
tempX = circular(unlist(plotX[[t]]))
tempY = circular(unlist(plotY[[t]]))
tempZ = circular(unlist(plotZ[[t]]))
resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
unifPhi_df <- as.data.frame(cbind(t=t, circleUniformPhi=circleUniformPhi))
plotX_df[[t]] <- as.data.frame(cbind(t=t, plotX=unlist(plotX[[t]])))
plotY_df[[t]] <- as.data.frame(cbind(t=t, plotY=unlist(plotY[[t]])))
plotZ_df[[t]] <- as.data.frame(cbind(t=t, plotZ=unlist(plotZ[[t]])))
pl_scene <- pl_scene %>% add_trace(data=unifPhi_df, showlegend=FALSE,
x = ~((t-epsilon)*cos(circleUniformPhi)),
y = ~((t-epsilon)*sin(circleUniformPhi)), z=0,
name=paste0("Time=",t), line=list(color='gray'),
mode = 'lines', opacity=0.3) %>%
add_markers(data=plotX_df[[t]], x=~(t*cos(plotX)), y=~(t*sin(plotX)), z=0,
type='scatter3d', name=paste0("X: t=",t),
marker=list(color='green'), showlegend=FALSE,
mode = 'markers', opacity=0.3) %>%
add_markers(data=plotY_df[[t]], x=~((t+epsilon)*cos(plotY)),
y=~((t+epsilon)*sin(plotY)), z=0-epsilon, showlegend=FALSE,
type='scatter3d', name=paste0("Y: t=",t),
marker=list(color='blue'),
mode = 'markers', opacity=0.3) %>%
add_markers(data=plotZ_df[[t]], x=~((t+2*epsilon)*cos(plotZ)),
y=~((t+2*epsilon)*sin(plotZ)), z=0+epsilon, showlegend=FALSE,
type='scatter3d', name=paste0("Z: t=",t),
marker=list(color='red'),
mode = 'markers', opacity=0.3)
}
means_df <- as.data.frame(cbind(t = c(1:length(r)),
plotX_means=unlist(lapply(plotX, mean)),
plotY_means=unlist(lapply(plotY, mean)),
plotZ_means=unlist(lapply(plotZ, mean))))
pl_scene <- pl_scene %>%
# add averaged (denoised) phase trajectories
add_trace(data=means_df, x=~(t*cos(plotX_means)),
y=~(t*sin(plotX_means)), z=0,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='green', width=15), opacity=0.8) %>%
add_trace(data=means_df, x=~(t*cos(plotY_means)),
y=~(t*sin(plotY_means)), z=0-epsilon,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='blue', width=15), opacity=0.8) %>%
add_trace(data=means_df, x=~(t*cos(plotZ_means)),
y=~(t*sin(plotZ_means)), z=0+epsilon,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='red', width=15), opacity=0.8) %>%
add_trace(x=0, y=0, z=c(-2,2), name="Space", showlegend=FALSE,
line=list(color='gray', width=15), opacity=0.8) %>%
layout(title="Pseudo Spacekime (1D Space, 2D Kime) Kime-Phase Sampling and Foliation",
scene = list(xaxis=list(title="Kappa1"), yaxis=list(title="Kappa2"),
zaxis=list(title="Space"))) %>% hide_colorbar()
pl_sceneThe intrinsic process error \(\varepsilon_m\sim N(0,\sigma^2V)\) reflects classical variability in the underlying process distribution. The more subtle error-tensor \(\varepsilon_g\sim N(0,R_g)\) depends on the variance-covariance matrix of the group-level errors \(R_g\) and can be interpreted as repeated-measures variability due to the kime-phase distribution.
Hence, the kime-phases can be used to model or track repeated-measure variations. Recall this kime-phase simulation we showed above demonstrating alternative kime-phase distributions coupled with corresponding mean-phase-aggregators (solid radial lines) representing real observable (deterministic, not stochastic) measurable quantities. The latter avoid the intrinsic kime-phase distribution noise due to random sampling through the phase-space distribution \(\Phi_{[-\pi,\pi)}^{(t)}\) when taking real measurements depending de facto on some phase-aggregating kernels that smooth the intractable (unobservable) and intrinsically-random kime-phases.
In applied sciences, sample statistics, such as the sample mean, variance, percentiles, range, interquartile range, etc., are always well-defined. Their values are not known at any given time until the actual (finite) sample \(\{x_i\}_{i=1}^n\) is collected (observations are recorded). For instance the sample mean statistic is the arithmetic average of all observations, \(\bar{x}_n=\frac{1}{n}\sum_{i-1}^n {x_i}\), yet, not knowing the value of the sample-mean prior to collecting the data is simply due to lack of evidence, since we have a closed form expression of how to obtain the sample average statistic estimating the population (physical system’s) expected mean response, \(\mu_X\ \underbrace{\leftarrow}_{{n\to\infty}}\ \bar{x}_n\).
Quantum systems interact in ways that can be explained with superpositions of different discrete base-states, and quantum system measurements yield statistical results corresponding to any one of the possible states appearing at random.
Similarly, analytical systems interact in ways that can be explained with superpositions of different, discrete, and finite samples. Correspondingly, analytical system measurements yield statistics corresponding to any one of the specific sample-statistic outcomes, which vary between samples and experiments. However, the first two fundamental laws of probability theory governing all statistical inference (CLT and LLN) yield that this between sample variability of sample-statistics is expected to decrease rapidly with the increases of the sample size. For instance, IID samples \(\{x_i\}_{i=1}^n\sim \mathcal{D}(\mu_X,\sigma_X^2)\), drawn from most nice distributions \(\mathcal{D}\) with mean \(\mu_X\) and variance \(\sigma_X^2\), yield sample-means \(\bar{x}_n\) whose sampling distribution converges (in distribution) to a Normal distribution, an asymptotic limiting distribution with rapidly decaying variance as the sample size increases
\[\underbrace{\bar{x}_n }_{Sample\\ average}\ \ \sim\ \underbrace{\bar{\mathcal{D}}(\overbrace{\mu_{\bar{x}_n}}^{mean}, \overbrace{\sigma^2_{\bar{x}_n}}^{variance}) }_{Sampling\\ distribution}\ \ \ \overbrace{\underbrace{\longrightarrow}_{n\to\infty}}^{Convergence\\ in\ distribution} \underbrace{\mathcal{N}\left(\mu_X,\frac{\sigma_X^2}{n}\right)}_{Asymptotic\\ distribution}\ .\]
Let’s try to explicate the duality between sampling from the kime-phase distribution and the complementary sampling of the process state space (spatiotemporal sampling).
Suppose the process \(X\sim F_X\), where \(F_X\) is an arbitrary distribution with a mean \(\mu=\mu_X\) and variance \(\sigma^2=\sigma^2_X\). If \(\{X_1,\cdots ,X_n\}\sim F_X\) is an IID sample the CLT suggests that under mild assumptions
\[\bar{X}\underset{n\to\infty}{\overset{d}{\longrightarrow}} N\left(\mu,\frac{\sigma^2}{n}\right )\ .\]
In repeated spatiotemporal sampling, \(\forall\ n\), we observe \(K\) (longitudinal) sample means \(\{\bar{X}_{n1},\bar{X}_{n2},\cdots, \bar{X}_{nK}\}\), where \(\forall\ 1\leq k_o\leq K\) each IID sample \(\{X_{1k_o},X_{2k_o},\cdots ,X_{nk_o}\}\sim F_X\) (across the kime-phase space) allows us to compute the \(k_o^{th}\) sample mean \(\bar{X}_{nk_o}\).
Observe that computing \(\bar{X}_{nk_o}\) from the \(k_o^{th}\) sample \(\{X_{1k_o},\cdots ,X_{nk_o}\}\) is identical to directly sampling \(\bar{X}_{nk_o}\) from the distribution \(F_{\bar{X}_n}\). Let’s reflect on this set up which clearly involves \(2\) independent sampling strategies.
Note that the distribution \(F_{\bar{X}_n}\) need not be Normal, it may be, but in general it won’t be exactly Normal. For instance suppose our spatiotemporal sampling scheme involves exponential distribution, i.e., \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim f_X\equiv {\text{Exp}}(\lambda)\equiv\Gamma(1,\lambda)\) and \(\sum_{i=1}^n{X_{ik_o}}\sim \Gamma\left(n,\lambda\right )\),
\[f_{{\text{Exp}}{(\lambda)}}(x)=\lambda e^{-\lambda x}\ \ ; \ \mathbb{E}(X_{{\text{Exp}}{(\lambda)}})=\frac{1}{\lambda} \ \ ; \ Var(X_{{\text{Exp}}{(\lambda)}})= \frac{1}{\lambda^2}\ ; \\ f_{\Gamma(shape=\alpha,scale=\lambda)}(x)=\frac{x^{\alpha-1}\lambda^{\alpha}}{\Gamma(\alpha)}e^{-\lambda x}\ \ ; \ \mathbb{E}(X_{\Gamma(\alpha,\lambda)})=\alpha\lambda \ \ ; \ Var(X_{\Gamma(\alpha,\lambda)})= \alpha\lambda^2\ \ .\]
Based on this sample, the sampling distribution of \(\bar{X}_{nk_o}\equiv \frac{1}{n}\sum_{i=1}^n{X_{ik_o}}\sim Exp\left (\lambda\right )\equiv\Gamma\left(n,n\lambda\right )\), since a sum of exponential IID \(Exp(\lambda)\) variables is \(\Gamma\)-distributed, \(\sum_{i=1}^n{X_{ik_o}}\sim \Gamma\left(n,\lambda\right )\) and scaling a random variable \(Y=cX\) yields \(F_Y(x)=\frac{1}{c}f_X(\frac{x}{c})\), in our case the constant multiplier \(c=\frac{1}{n}\), and \(Y=\frac{1}{n}\sum_{i=1}^n{X_i}\).
Of course, as \(n\to\infty\), \(\Gamma\left(n,\frac{\gamma}{n}\right )\to N(0,1)\), yet for any fixed \(n\), the (gamma) \(\Gamma\) distribution is similar, but not identical to normal.
Prior to data acquisition, \(\bar{X}_{nk_o}\) is a random variable, once the observed data values are plugged in, it’s a constant. Hence, the sample mean random variable (a sum of \(n\) \(\sim F_X\) variables) \(\bar{X}_{nk_o}\) based on \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim F_X\), and the (single!) random variable \(Y\sim F_{\bar{X}_{nk_o}}\) represent exactly the same random variable.
Corollary: Random drawing of a second-order tensor \(X_{n\times K}\) (a random matrix) representing \(K\) (spatiotemporal or longitudinal) samples of \(n\) IID observations \(\{X_{1k_o},\cdots , X_{nk_o}\}_{k_o}\sim F_X\) and then computing the \(K\) corresponding means \(\bar{X}_{nk_o}=\frac{1}{n}\sum_{i=1}^n{X_{ik_o}}\) is equivalent to randomly drawing \(K\) samples (a random vector) directly from the sampling distribution of the mean \(F_{\bar{X}_{nk_o}}\).
Below is an example of a 3D Brownian motion/Wiener process as an example of a stochastic random walk process. In this example, the Wiener process is intentionally disturbed by random Poisson noise, which leads to occasional abrupt disruptions of the already stochastic process.
library(plotly)
# 3D Wiener Process
N=500 # Number of random walk steps
# Define the X, Y, and Z, coordinate displacements independently
# xdis = rnorm(N, 0 , 1)
# ydis = rnorm(N, 0 , 2)
# zdis = rnorm(N, 0 , 3)
# xdis = cumsum(xdis)
# ydis = cumsum(ydis)
# zdis = cumsum(zdis)
# To use simulated 3D MVN Distribution
dis <- mixtools::rmvnorm(N, mu=c(0,0,0),
sigma=matrix(c(1,0,0, 0,2,0, 0,0,3), ncol=3))
# aggregate the displacements to get the actual 3D Cartesian Coordinates
xdis = cumsum(dis[,1])
ydis = cumsum(dis[,2])
zdis = cumsum(dis[,3])
# add Poisson noise
at = rpois(N, 0.1)
for(i in c(1:N)) {
if(at[i] != 0) {
xdis[i] = xdis[i]*at[i]
ydis[i] = ydis[i]*at[i]
zdis[i] = ydis[i]*at[i]
}
}
# plot(xdis, ydis, type="l",
# main ="Brownian Motion in Two Dimension with Poisson Arrival Process",
# xlab="x displacement", ylab = "y displacement")
plot_ly(x=xdis, y=ydis, z=zdis, type="scatter3d", mode="markers+lines",
text=~c(1:N), hoverinfo='text',
marker=list(color='gray'), showlegend=F) %>%
# Emphasize the starting and ending points
add_markers(x=xdis[1], y=ydis[1], z=zdis[1], marker=list(size=20,color="green"),
text=paste0("Starting Node 1")) %>%
add_markers(x=xdis[N], y=ydis[N], z=zdis[N],
marker=list(size=20,color="red"), text=paste0("Ending Node ", N)) %>%
layout(title="3D Brownian Motion Simulation (MVN & Poisson)",
scene=list(xaxis=list(title="X displacement"),
yaxis=list(title="Y displacement"),
zaxis=list(title="Z displacement")))# 1D Wiener Process
# dis = rnorm(N, 0, 1);
# at = rpois(N,1)
# for(i in 1:N) {
# if(at[i] != 0){
# dis[i]= dis[i]*at[i]
# }
# }
# dis = cumsum(dis)
# plot(dis, type= "l",
# main= "Brownian Motion in One Dimension with Poisson Arrival Process",
# xlab="time", ylab="displacement")
# ub = 20; lb = -20
# xdis = rnorm(N, 0 ,1)
# xdis1 = rep(1,N)
# xdis1[1] = xdis[1]
# for(i in c(1:(N-1))){
# if(xdis1[i] + xdis[i+1] > ub) { xdis1[i+1] <- ub }
# else if(xdis1[i] + xdis[i+1] < lb) { xdis[i+1] = lb }
# else { xdis1[i+1] = xdis1[i] + xdis[i+1] }
# }
#
# plot(xdis1, type="l",main="Brownian Motion with bound in 1-dim", xlab="displacement",ylab="time")
# Compute the row Euclidean differences
df <- data.frame(cbind(xdis, ydis, zdis))
rowEuclidDistance <- dist(df)
plot_ly(z=as.matrix(rowEuclidDistance), type="heatmap") %>%
layout(title="Heatmap of Euclidean Distances between Consecutive Steps of Wiener Process")