fit_regression {nprotreg}

Fits a 3D Spherical Regression.

Description

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

Usage

fit_regression(
  evaluation_points,
  explanatory_points,
  response_points,
  concentration,
  weights_generator = weight_explanatory_points,
  number_of_expansion_terms = 1,
  number_of_iterations = 1,
  allow_reflections = FALSE
)

Arguments

evaluation_points

An n-by-3 matrix whose rows contain the Cartesian coordinates of the points at which the regression will be estimated.

explanatory_points

An m-by-3 matrix whose rows contain the Cartesian coordinates of the explanatory points used to calculate the regression estimators.

response_points

An m-by-3 matrix whose rows contain the Cartesian coordinates of the response points corresponding to the explanatory points.

concentration

A non negative scalar whose reciprocal value is proportional to the bandwidth applied while estimating a spherical regression model.

weights_generator

A function that, given a matrix of n evaluation points, returns an m-by-n matrix whose j-th column contains the weights assigned to the explanatory points while analyzing the j-th evaluation point. Defaults to weight_explanatory_points.

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 1 or 2. Defaults to 1.

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 1.

allow_reflections

A logical scalar value. If set to TRUE signals that reflections are allowed. Defaults to FALSE. It is ignored if number_of_expansion_terms is 2.

Details

Function weights_generator must be prototyped as having the following three arguments:

evaluation_points

a matrix whose n rows are the Cartesian coordinates of given evaluation points.

explanatory_points

a matrix whose m rows are the Cartesian coordinates of given explanatory points.

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 weights_generator will return a non NULL numerical m-by-n matrix whose j-th column contains the weights assigned to the explanatory points while analyzing the j-th evaluation point.

Function fit_regression supports parallel execution. To setup parallelization, you can exploit the doParallel package. Otherwise, fit_regression will be executed sequentially and, when called the first time, you will receive the following

## Warning: executing %dopar% sequentially: no parallel backend registered

This is completely safe and by design.

Value

A number_of_iterations-length vector of lists, with the s-th list having two components, fitted_response_points, an n-by-3 matrix whose rows contain the Cartesian coordinates of the fitted points at iteration s, and explanatory_points, an m-by-3 matrix whose rows contain the Cartesian coordinates of the points exploited as explanatory at iteration s.

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>.

See Also

Other Regression functions: cross_validate_concentration(), get_equally_spaced_points(), get_skew_symmetric_matrix(), simulate_regression(), simulate_rigid_regression(), weight_explanatory_points()

Examples

library(nprotreg)

# Create 100 equally spaced design points on the sphere.

number_of_explanatory_points <- 100

explanatory_points <- get_equally_spaced_points(
  number_of_explanatory_points
)

# Define the regression model, where the rotation for a given "point"
# is obtained from the exponential of a skew-symmetric matrix with the
# following components.

local_rotation_composer <- function(point) {
  independent_components <- (1 / 8) *
    c(exp(2.0 * point[3]), - exp(2.0 * point[2]), exp(2.0 * point[1]))
}

# Define an error term given by a small rotation, similarly defined
# from a skew-symmetric matrix with random entries.

local_error_sampler <- function(point) {
  rnorm(3, sd = .01)
}

# Generate the matrix of responses, using the regression model
# and the error model.

response_points <- simulate_regression(
  explanatory_points,
  local_rotation_composer,
  local_error_sampler
)

# Create some "test data" for which the response will be predicted.

evaluation_points <- rbind(
  cbind(.5, 0, .8660254),
  cbind(-.5, 0, .8660254),
  cbind(1, 0, 0),
  cbind(0, 1, 0),
  cbind(-1, 0, 0),
  cbind(0, -1, 0),
  cbind(.5, 0, -.8660254),
  cbind(-.5, 0, -.8660254)
)

# Define a weight function for nonparametric fit.

weights_generator <- weight_explanatory_points

# Set the concentration parameter.

concentration <- 5

# Or obtain this by cross-validation: see
# the `cross_validate_concentration` function.

# Fit regression.

fitted_model <- fit_regression(
  evaluation_points,
  explanatory_points,
  response_points,
  concentration,
  weights_generator,
  number_of_expansion_terms = 1,
  number_of_iterations = 2
)

# Extract the point corresponding to the
# second evaluation point fitted at
# the first iteration.

cat("Point fitted at iteration 1 corresponding to the second evaluation point: \n")
cat(fitted_model[[1]]$fitted_response_points[2, ], "\n")

## Not run: 
# Create some plots to view the results.

# 3D plot.

library(rgl)

plot3d(
  explanatory_points,
  type = "n",
  xlab = "x",
  ylab = "y",
  zlab = "z",
  box = TRUE,
  axes = TRUE
)
spheres3d(0, 0, 0, radius = 1, lit = FALSE, color = "white")
spheres3d(0, 0, 0, radius = 1.01, lit = FALSE, color = "black", front = "culled")
text3d(c(0, 0, 1), text = "N", adj = 0)

ll <- 10
vv1 <- (ll - (0:(ll))) / ll
vv2 <- 1 - vv1
plot3d(explanatory_points, add = TRUE, col = 2)
for (i in 1:dim(explanatory_points)[1]) {
  m <- outer(vv1, explanatory_points[i,], "*") +
    outer(vv2, response_points[i,], "*")
  m <- m / sqrt(apply(m ^ 2, 1, sum))
  lines3d(m, col = 3)
}

plot3d(evaluation_points, add = TRUE, col = 4)

for (i in 1:dim(evaluation_points)[1]) {
  m <- outer(vv1, evaluation_points[i,], "*") +
    outer(vv2, fitted_model[[1]]$fitted_response_points[i,], "*")
  m <- m / sqrt(apply(m ^ 2, 1, sum))
  lines3d(m, col = 1)
}

# 2D plot.

explanatory_spherical_coords <- convert_cartesian_to_spherical(explanatory_points)
response_spherical_coords <- convert_cartesian_to_spherical(response_points)

plot(
  x = explanatory_spherical_coords[, 1],
  y = explanatory_spherical_coords[, 2],
  pch = 20,
  cex = .7,
  col = 2,
  xlab = "longitude",
  ylab = "latitude"
)

for (i in 1:dim(explanatory_spherical_coords)[1]) {
  column <- 1
  if ((explanatory_spherical_coords[i, 1] - response_spherical_coords[i, 1]) ^ 2 +
      (explanatory_spherical_coords[i, 2] - response_spherical_coords[i, 2]) ^ 2 > 4)
        column <- "grey"
  lines(
    c(explanatory_spherical_coords[i, 1], response_spherical_coords[i, 1]),
    c(explanatory_spherical_coords[i, 2], response_spherical_coords[i, 2]),
    col = column
  )
}

evaluation_spherical_coords <- convert_cartesian_to_spherical(
  evaluation_points
)

fitted_response_spherical_coords <- convert_cartesian_to_spherical(
  fitted_model[[1]]$fitted_response_points
)

points(
  x = evaluation_spherical_coords[, 1],
  y = evaluation_spherical_coords[, 2],
  pch = 20,
  cex = .7,
  col = 4
)

for (i in 1:dim(evaluation_spherical_coords)[1]) {
  column <- 3
  if ((evaluation_spherical_coords[i, 1] - fitted_response_spherical_coords[i, 1]) ^ 2 +
      (evaluation_spherical_coords[i, 2] - fitted_response_spherical_coords[i, 2]) ^ 2 > 4)
        column <- "grey"
  lines(
    c(evaluation_spherical_coords[i, 1], fitted_response_spherical_coords[i, 1]),
    c(evaluation_spherical_coords[i, 2], fitted_response_spherical_coords[i, 2]),
    col = column
  )
}


## End(Not run)