fit_regression: Fits a 3D Spherical Regression.
In nprotreg: Nonparametric Rotations for Sphere-Sphere Regression

fit_regression

R Documentation

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

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)

nprotreg documentation built on Sept. 28, 2023, 9:06 a.m.