Nothing
#' Fits a 3D Spherical Regression.
#'
#' Returns 3D spherical points obtained by locally rotating
#' the specified evaluation
#' points, given an approximated model for local rotations
#' and a weighting scheme for the observed data set.
#' This function implements the method for sphere-sphere regression
#' proposed by Di Marzio et al. (2018).
#'
#' Function \code{weights_generator} must be prototyped as having the
#' following three arguments:
#' \describe{
#' \item{\code{evaluation_points}}{a matrix whose \emph{n} rows are the Cartesian coordinates of
#' given evaluation points.}
#' \item{\code{explanatory_points}}{a matrix whose \emph{m} rows are the Cartesian coordinates of
#' given explanatory points.}
#' \item{\code{concentration}}{A non negative scalar whose reciprocal value
#' is proportional to the bandwidth applied while estimating
#' a spherical regression model.}
#' }
#' It is also expected that \code{weights_generator} will return
#' a non \code{NULL} numerical \emph{m}-by-\emph{n} matrix whose \emph{j}-th column contains
#' the weights assigned to the explanatory points while analyzing the
#' \emph{j}-th evaluation point.
#'
#' Function \code{fit_regression} supports parallel execution.
#' To setup parallelization, you can exploit the
#' \href{https://cran.r-project.org/package=doParallel}{doParallel} package.
#' Otherwise, \code{fit_regression} will be executed sequentially and, when called the
#' first time, you will receive the following
#'
#' \code{## Warning: executing \%dopar\% sequentially: no parallel backend registered}
#'
#' This is completely safe and by design.
#'
#' @param evaluation_points An \emph{n}-by-3 matrix whose rows contain
#' the Cartesian coordinates of the points at which the regression
#' will be estimated.
#' @param explanatory_points An \emph{m}-by-3 matrix whose rows contain
#' the Cartesian coordinates of the explanatory points used to
#' calculate the regression estimators.
#' @param response_points An \emph{m}-by-\emph{3} matrix whose rows contain
#' the Cartesian coordinates of the response points corresponding
#' to the explanatory points.
#' @param concentration A non negative scalar whose reciprocal value
#' is proportional to the bandwidth applied while estimating
#' a spherical regression model.
#' @param weights_generator A function that, given a matrix of \emph{n}
#' evaluation points, returns an \emph{m}-by-\emph{n} matrix whose
#' \emph{j}-th column contains
#' the weights assigned to the explanatory points while analyzing the
#' \emph{j}-th evaluation point. Defaults to \code{\link{weight_explanatory_points}}.
#' @param number_of_expansion_terms The number of terms to be included
#' in the expansion of the matrix exponential applied while
#' approximating a local rotation matrix. Must be \code{1} or \code{2}.
#' Defaults to \code{1}.
#' @param number_of_iterations The number of
#' rotation fitting steps to be executed.
#' At each step, the points estimated during the previous step
#' are exploited as the current explanatory points. Defaults to \code{1}.
#' @param allow_reflections A logical scalar value. If set to \code{TRUE} signals
#' that reflections are allowed. Defaults to \code{FALSE}. It is ignored if
#' \code{number_of_expansion_terms} is \code{2}.
#' @return A \code{number_of_iterations}-length vector of lists, with the \code{s}-th
#' list having two components,
#' \code{fitted_response_points}, an \emph{n}-by-\emph{3} matrix whose rows contain
#' the Cartesian coordinates of the fitted points at iteration \code{s}, and
#' \code{explanatory_points}, an \emph{m}-by-\emph{3} matrix whose rows contain
#' the Cartesian coordinates of the points exploited as explanatory at
#' iteration \code{s}.
#' @export
#' @family Regression functions
#' @example examples/example_fit_regression.R
#' @references Marco Di Marzio, Agnese Panzera & Charles C. Taylor (2018)
#' Nonparametric rotations for sphere-sphere regression,
#' Journal of the American Statistical Association,
#' <doi:10.1080/01621459.2017.1421542>.
fit_regression <- function(
evaluation_points,
explanatory_points,
response_points,
concentration,
weights_generator = weight_explanatory_points,
number_of_expansion_terms = 1,
number_of_iterations = 1,
allow_reflections = FALSE
) {
# Input validation ---
parameter_validators$evaluation_points$validate(evaluation_points, -4)
parameter_validators$explanatory_points$validate(explanatory_points, -4)
parameter_validators$response_points$validate(response_points, -4)
parameter_validators$response_points$check_shared_number_of_rows(
response_points,
explanatory_points,
parameter_validators$explanatory_points$name,
-3
)
parameter_validators$concentration$validate(concentration, -4)
fnc_typed_parameter_validators$weights_generator$validate(weights_generator, -4)
parameter_validators$number_of_expansion_terms$validate(number_of_expansion_terms, -4)
parameter_validators$number_of_iterations$validate(
number_of_iterations,
-4
)
parameter_validators$allow_reflections$validate(allow_reflections, -4)
# Implementation ---
local_rotation_modeler <- switch(
number_of_expansion_terms,
one_term_local_rotation_modeler,
two_term_local_rotation_modeler
)
number_of_evaluation_points <- nrow(evaluation_points)
transposed_response_points <- t(response_points)
explanatory_point_weights <- weights_generator(
evaluation_points,
explanatory_points,
concentration
)
explanatory_points_sequence <-
get_explanatory_points_sequence_par(
explanatory_points,
transposed_response_points,
concentration,
weights_generator,
number_of_iterations,
local_rotation_modeler,
allow_reflections
)
fit_info <- vector(mode = "list", length = number_of_iterations)
for (s in 1:number_of_iterations) {
fit_info[[s]] <- list(
fitted_response_points = matrix(
nrow = number_of_evaluation_points,
ncol = 3),
explanatory_points = explanatory_points_sequence[[s]]$explanatory_points
)
}
for (i in 1:number_of_evaluation_points) {
iterative_fitted_response_points <- fit_iterative_regression(
evaluation_points[i, ],
explanatory_points_sequence,
transposed_response_points,
explanatory_point_weights[, i],
local_rotation_modeler,
allow_reflections
)
for (s in 1:number_of_iterations) {
fit_info[[s]]$fitted_response_points[i, ] <- iterative_fitted_response_points[s, ]
}
}
fit_info
#----
}
# Implements and generalizes Algorithm 1 for iterative rotation fitting.
#
# This is a generalization of the algorithm
# Di Marzio et al. (2018) refer to as 1.
# Here it is parameterized w.r.t.
# a generic local rotation modeler.
#
# This function expects the response points in a transposed form,
# in order to enhance the performance of rotation estimations
# for one-term models.
#
# @param evaluation_point A row containing
# the Cartesian coordinates of the point at which the regression
# will be estimated.
# @param explanatory_points_sequence A \code{number_of_iterations}-length
# vector of lists having the single component \code{explanatory_points}.
# In the list at position \code{s}, it is the \emph{m}-by-3 matrix whose rows contain
# the Cartesian coordinates of the explanatory points used to
# calculate the regression estimators at iteration \code{s}.
# @param transposed_response_points A \emph{3}-by-\emph{m} matrix whose columns
# contain the Cartesian coordinates of the response points corresponding
# to the explanatory points.
# @param evaluation_specific_explanatory_point_weights
# An \eqn{m}-by-\eqn{1} vector whose entries
# store the weights assigned to the observed
# \code{explanatory_points} while estimating the regression
# at the given \code{evaluation_point}.
# @param local_rotation_modeler A model for a local rotation.
# It must be prototyped as having three arguments, \code{evaluation_point},
# \code{explanatory_points}, and \code{transposed_response_points},
# and returning a 3D rotation matrix.
# @param allow_reflections A logical scalar value. If set to \code{TRUE} signals
# that reflections are allowed. Defaults to \code{FALSE}. It is ignored if
# \code{number_of_expansion_terms} is \code{2}.
#
# @return A \code{number_of_iterations}-by-\eqn{3} matrix whose \code{s}-th row
# contains the response obtained applying the given rotation model
# to the specified \code{evaluation_point} at iteration \code{s}.
fit_iterative_regression <- function(
evaluation_point,
explanatory_points_sequence,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
local_rotation_modeler,
allow_reflections
) {
number_of_iterations <- length(explanatory_points_sequence)
iterative_fitted_response_points <- matrix(
nrow = number_of_iterations,
ncol = 3
)
cumulative_rotation <- diag(3)
for (s in 1:number_of_iterations) {
additional_rotation <- local_rotation_modeler(
evaluation_point,
explanatory_points_sequence[[s]]$explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
allow_reflections)
cumulative_rotation <- additional_rotation %*% cumulative_rotation
iterative_fitted_response_points[s, ] <- cumulative_rotation %*% evaluation_point
}
iterative_fitted_response_points
}
# Gets the sequence of explanatory points matrices needed to
# to implement Algorithm 1 for iterative rotation fitting.
#
# This function supports parallel execution (parallel version).
#
# This function expects the response points in a transposed form,
# in order to enhance the performance of rotation estimations
# for one-term models.
#
# @param explanatory_points An \emph{m}-by-3 matrix whose rows contain
# the Cartesian coordinates of the explanatory points used to
# calculate the regression estimators.
# @param transposed_response_points A \emph{3}-by-\emph{m} matrix whose columns
# contain the Cartesian coordinates of the response points corresponding
# to the explanatory points.
# @param concentration A non negative scalar whose reciprocal value
# is proportional to the bandwidth applied while estimating
# a spherical regression model.
# @param weights_generator A function that, given a matrix of \emph{n}
# evaluation points, returns an \emph{m}-by-\emph{n} matrix whose
# \emph{j}-th column contains
# the weights assigned to the explanatory points while analyzing the
# \emph{j}-th evaluation point. Defaults to \code{\link{weight_explanatory_points}}.
# @param number_of_iterations The number of
# rotation fitting steps to be executed.
# At each step, the points estimated during the previous step
# are exploited as the current explanatory points. Defaults to \code{1}.
# @param local_rotation_modeler A model for a local rotation.
# It must be prototyped as having three arguments, \code{evaluation_point},
# \code{explanatory_points}, and \code{transposed_response_points},
# and returning a 3D rotation matrix.
# @param allow_reflections A logical scalar value. If set to \code{TRUE} signals
# that reflections are allowed. Defaults to \code{FALSE}. It is ignored if
# \code{number_of_expansion_terms} is \code{2}.
#
# @return A \code{number_of_iterations}-length vector of lists, with the \code{s}-th
# list having a single component,
# \code{explanatory_points}, an \emph{m}-by-\emph{3} matrix whose rows contain
# the Cartesian coordinates of the points exploited as explanatory at
# iteration \code{s}.
get_explanatory_points_sequence_par <- function(
explanatory_points,
transposed_response_points,
concentration,
weights_generator,
number_of_iterations,
local_rotation_modeler,
allow_reflections
) {
explanatory_points_sequence <- vector(
mode = "list",
length = number_of_iterations
)
explanatory_points_sequence[[1]]$explanatory_points <- explanatory_points
if (number_of_iterations > 1) {
evaluation_specific_explanatory_point_weights <-
weights_generator(
explanatory_points,
explanatory_points,
concentration
)
number_of_explanatory_points <- nrow(explanatory_points)
for (s in 2:number_of_iterations) {
iteration_evaluation_points <- explanatory_points_sequence[[s - 1]]$explanatory_points
i <- NULL
iteration_explanatory_points <-
foreach(i = 1:number_of_explanatory_points, .combine = "rbind") %dopar% {
iteration_fitted_response_point <- fit_iterative_regression(
iteration_evaluation_points[i, ],
explanatory_points_sequence[s - 1],
transposed_response_points,
evaluation_specific_explanatory_point_weights[, i],
local_rotation_modeler,
allow_reflections)
iteration_fitted_response_point
}
explanatory_points_sequence[[s]]$explanatory_points <- iteration_explanatory_points
}
}
explanatory_points_sequence
}
# Models Rotations With One-Term Local Expansions.
#
# This is an implementation of the one-term local rotation
# model proposed by Di Marzio et al. (2018).
#
# @param evaluation_point A row containing
# the Cartesian coordinates of the point at which the regression
# will be estimated.
# @param explanatory_points An \emph{m}-by-3 matrix whose rows contain
# the Cartesian coordinates of the explanatory points used to
# calculate the regression estimators.
# @param transposed_response_points A \emph{3}-by-\emph{m} matrix whose columns
# contain the Cartesian coordinates of the response points corresponding
# to the explanatory points.
# @param evaluation_specific_explanatory_point_weights
# An \eqn{m}-by-\eqn{1} vector whose entries
# store the weights assigned to the observed
# \code{explanatory_points} while estimating the regression
# at the given \code{evaluation_point}.
# @param allow_reflections A logical scalar value. If set to \code{TRUE} signals
# that reflections are allowed. Defaults to \code{FALSE}. It is ignored if
# \code{number_of_expansion_terms} is \code{2}.
#
# @return The modeled rotation matrix local to the specified
# \code{evaluation_point}.
one_term_local_rotation_modeler <- function(
evaluation_point,
explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
allow_reflections
) {
svd_info <- svd(
transposed_response_points %*%
diag(evaluation_specific_explanatory_point_weights) %*%
explanatory_points
)
u <- svd_info$u
transposed_v <- t(svd_info$v)
local_rotation <- u %*% transposed_v
if (!allow_reflections) {
local_rotation_det <- det(local_rotation)
local_rotation <- u %*% diag(c(1, 1, local_rotation_det)) %*% transposed_v
}
local_rotation
}
# Models Rotations With Two-Term Local Expansions.
#
# This is an implementation of the two-term local rotation
# model proposed by Di Marzio et al. (2018).
#
# @param evaluation_point A row containing
# the Cartesian coordinates of the point at which the regression
# will be estimated.
# @param explanatory_points An \emph{m}-by-3 matrix whose rows contain
# the Cartesian coordinates of the explanatory points used to
# calculate the regression estimators.
# @param transposed_response_points A \emph{3}-by-\emph{m} matrix whose columns
# contain the Cartesian coordinates of the response points corresponding
# to the explanatory points.
# @param evaluation_specific_explanatory_point_weights
# An \eqn{m}-by-\eqn{1} vector whose entries
# store the weights assigned to the observed
# \code{explanatory_points} while estimating the regression
# at the given \code{evaluation_point}.
# @param allow_reflections A logical scalar value. If set to \code{TRUE} signals
# that reflections are allowed. Defaults to \code{FALSE}. It is ignored if
# \code{number_of_expansion_terms} is \code{2}.
#
# @return The modeled rotation matrix local to the specified
# \code{evaluation_point}.
two_term_local_rotation_modeler <- function(
evaluation_point,
explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
allow_reflections
) {
# Function crude2 ---
# Computes The 3D Weighted Least Squares As A Function Of 12 components.
#
# @param components A 12-length vector interpreted as the argument
# of a function whose value is
# the Weighted Least Squares expression.
#
# @return The WLS for the specified \code{components}.
#
# @inheritParams fit_iterative_regression
final_objective_function <- function(
evaluation_point,
explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
components
) {
# This is equivalent to function 'crude2' in legacy code.
r <- components # matrix(components, 3, 4)
dim(r) <- c(3, 4)
w <- evaluation_specific_explanatory_point_weights
x <- explanatory_points
y <- t(transposed_response_points)
transposed_x <- t(x)
number_of_points <- nrow(x)
id <- diag(3)
yhat <- matrix(0.0, nrow = number_of_points, ncol = 3)
dd <- rbind(1, evaluation_point - transposed_x)
r1 <- r %*% dd
for (i in 1:number_of_points) {
current_x <- transposed_x[, i]
yhat[i, ] <- current_x
if (!all(r1[, i] == 0)) {
s <- get_skew_symmetric_matrix(r1[, i])
s_norm <- sqrt(s[1, 2] ^ 2 + s[1, 3] ^ 2 + s[2, 3] ^ 2)
s <- s / s_norm
local_rotation <- id + sin(s_norm) * s + (1 - cos(s_norm)) * s %*% s
yhat[i, ] <- local_rotation %*% current_x
}
}
value <- sum(w * (y - yhat) ^ 2)
}
# Implementation ---
r <- one_term_local_rotation_modeler(
evaluation_point,
explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights,
allow_reflections
)
s0 <- logm(r)
rr <- c(-s0[2, 3], s0[1, 3], -s0[1, 2])
pp <- c(rr, 0, 0, 0, 0, 0, 0, 0, 0, 0)
nn1 <- nlm(
final_objective_function,
p = pp,
evaluation_point = evaluation_point,
explanatory_points = explanatory_points,
transposed_response_points = transposed_response_points,
evaluation_specific_explanatory_point_weights =
evaluation_specific_explanatory_point_weights
)
s <- get_skew_symmetric_matrix(nn1$est[1:3])
s_norm <- sqrt(s[1, 2] ^ 2 + s[1, 3] ^ 2 + s[2, 3] ^ 2)
s <- s / s_norm
local_rotation <- diag(3) + sin(s_norm) * s + (1 - cos(s_norm)) * s %*% s
# ----
}
# Computes Weighted Least Squares For 3D Regressions.
#
# This is the objective function of the problem
# described in Equation (9).
#
# @param local_rotation A 3-by-3 rotation matrix.
# @param explanatory_points An \emph{m}-by-3 matrix whose rows contain
# the Cartesian coordinates of the explanatory points used to
# calculate the regression estimators.
# @param transposed_response_points A \emph{3}-by-\emph{m} matrix whose columns
# contain the Cartesian coordinates of the response points corresponding
# to the explanatory points.
# @param evaluation_specific_explanatory_point_weights
# An \eqn{m}-by-\eqn{1} vector whose entries
# store the weights assigned to the observed
# \code{explanatory_points} while estimating the regression
# at the given \code{evaluation_point}.
# @return The weighted least squares corresponding
# to the specified parameters.
weighted_least_squares <- function(
local_rotation,
explanatory_points,
transposed_response_points,
evaluation_specific_explanatory_point_weights
) {
r <- local_rotation
w <- evaluation_specific_explanatory_point_weights
x <- explanatory_points
yT <- transposed_response_points
hT <- yT - t(x %*% r)
sum(diag(hT %*% diag(w) %*% t(hT)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.