# R/special_sphere.R In Riemann: Learning with Data on Riemannian Manifolds

#### Documented in sphere.geo2xyzsphere.runifsphere.utestsphere.xyz2geo

## SOME SPECIAL FUNCTIONS ON SPHERE
#  (01) sphere.runif
#  (02) sphere.utest
#  (03) sphere.convert

# (01) sphere.runif ============================================================
#' Generate Uniform Samples on Sphere
#'
#' It generates \eqn{n} random samples from \eqn{\mathcal{S}^{p-1}}. For convenient
#' usage of users, we provide a number of options in terms of the return type.
#'
#' @param n number of samples to be generated.
#' @param p original dimension (of the ambient space).
#' @param type return type; \describe{
#' \item{\code{"list"}}{a length-\eqn{n} list of length-\eqn{p} vectors.}
#' \item{\code{"matrix"}}{a \eqn{(n\times p)} where rows are unit vectors.}
#' \item{\code{"riemdata"}}{a S3 object. See \code{\link{wrap.sphere}} for more details (\emph{Default}).}
#' }
#'
#' @return an object from one of the above by \code{type} option.
#'
#' @examples
#' #-------------------------------------------------------------------
#' #                       Draw Samples on Sphere
#' #
#' # Multiple return types on S^4 in R^5
#' #-------------------------------------------------------------------
#' dat.list = sphere.runif(n=10, p=5, type="list")
#' dat.matx = sphere.runif(n=10, p=5, type="matrix")
#' dat.riem = sphere.runif(n=10, p=5, type="riemdata")
#'
#' @references
#' \insertRef{chikuse_statistics_2003}{Riemann}
#'
#' @concept sphere
#' @export
sphere.runif <- function(n, p, type=c("list","matrix","riemdata")){
# PREPROCESSING
# parameters
N = round(n)
p = round(p) # S^{p-1} in R^p

# return object type
if (missing(type)){
retype = "riemdata"
} else {
retype = match.arg(tolower(type), c("list","matrix","riemdata"))
}

# GENERATION
outmat = runif_sphere(N,p)
if (all(retype=="matrix")){
return(outmat)
} else {
outobj = wrap.sphere(outmat)
if (all(retype=="list")){
return(outobj$data) } else { return(outobj) } } } # (02) sphere.utest ============================================================ #' Test of Uniformity on Sphere #' #' Given \eqn{N} observations \eqn{\lbrace X_1, X_2, \ldots, X_M \brace} on #' \eqn{\mathcal{S}^{p-1}}, it tests whether the data is distributed uniformly #' on the sphere. #' #' @param spobj a S3 \code{"riemdata"} class for \eqn{N} Sphere-valued data. #' @param method (case-insensitive) name of the test method containing \describe{ #' \item{\code{"Rayleigh"}}{original Rayleigh statistic.} #' \item{\code{"RayleighM"}}{modified Rayleigh statistic with better order of error.} #' } #' #' @return a (list) object of \code{S3} class \code{htest} containing: \describe{ #' \item{statistic}{a test statistic.} #' \item{p.value}{\eqn{p}-value under \eqn{H_0}.} #' \item{alternative}{alternative hypothesis.} #' \item{method}{name of the test.} #' \item{data.name}{name(s) of provided sample data.} #' } #' #' #' @examples #' #------------------------------------------------------------------- #' # Compare Rayleigh's original and modified versions of the test #' #------------------------------------------------------------------- #' # Data Generation #' myobj = sphere.runif(n=100, p=5, type="riemdata") #' #' # Compare 2 versions : Original vs Modified Rayleigh #' sphere.utest(myobj, method="rayleigh") #' sphere.utest(myobj, method="rayleighm") #' #' #' @references #' \insertRef{chikuse_statistics_2003}{Riemann} #' #' \insertRef{mardia_directional_1999}{Riemann} #' #' @seealso \code{\link{wrap.sphere}} #' @concept sphere #' @export sphere.utest <- function(spobj, method=c("Rayleigh","RayleighM")){ ## CHECK INPUT check_inputmfd(spobj, "sphere.utest") DNAME = deparse(substitute(spobj)) # borrowed from HDtest method.now = ifelse(missing(method),"rayleigh", match.arg(tolower(method), c("rayleigh","rayleighm"))) ## COMPUTE # prepare data in 3d array data3d <- aux_rvec2array3d(spobj) output <- switch(method.now, rayleigh = sp.utest.Rayleigh(data3d, DNAME, is.original = TRUE), rayleighm = sp.utest.Rayleigh(data3d, DNAME, is.original = FALSE)) return(output) } #' @keywords internal #' @noRd sp.utest.Rayleigh <- function(x, dname, is.original=TRUE){ # Take the version from RiemStiefel {will be deprecated} p = dim(x)[1] # r-frame in R^p with n observations r = dim(x)[2] n = dim(x)[3] xbar = array(0,c(p,r)) for (i in 1:p){ for (j in 1:r){ xbar[i,j] = base::mean(as.vector(x[i,j,])) } } S = p*n*sum(diag(t(xbar)%*%xbar)) if (is.original){ hname = "Rayleigh Test of Uniformity on Sphere" thestat = S pvalue = stats::pchisq(thestat, df=as.integer(p*r), lower.tail=FALSE) } else { hname = "Modified Rayleigh Test of Uniformity on Sphere" term1 = 1/(2*n) term2 = 1 - (S/((p*r) + 2)) thestat = S*(1 - term1*term2) pvalue = stats::pchisq(thestat, df=as.integer(p*r), lower.tail=FALSE) } # COMPUTATION : DETERMINATION Ha = paste("data is not uniformly distributed on ",p-1,"-sphere.",sep="") names(thestat) = "statistic" res = list(statistic=thestat, p.value=pvalue, alternative = Ha, method=hname, data.name = dname) class(res) = "htest" return(res) } # (03) sphere.convert ========================================================= # https://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates #' Convert between Cartesian Coordinates and Geographic Coordinates #' #' In geospatial data analysis, it is common to consider locations on the Earth as #' data. These locations, usually provided by latitude and longitude, are not directly #' applicable for spherical data analysis. We provide two functions - \code{sphere.geo2xyz} and \code{sphere.xyz2geo} - #' that convert geographic coordinates in longitude/latitude into a unit-norm vector on \eqn{\mathcal{S}^2}, and vice versa. #' As a convention, latitude and longitude are represented as \emph{decimal degrees}. #' #' @param lat latitude (in decimal degrees). #' @param lon longitude (in decimal degrees). #' @param xyz a unit-norm vector in \eqn{\mathcal{S}^{2}}. #' #' @return transformed data. #' #' @examples #' ## EXAMPLE DATA WITH POPULATED US CITIES #' data(cities) #' #' ## SELECT ALBUQUERQUE #' geo = cities$coord[1,]
#' xyz = cities\$cartesian[1,]
#'
#' ## CHECK TWO INPUT TYPES AND THEIR CONVERSIONS
#' sphere.geo2xyz(geo[1], geo[2])
#' sphere.xyz2geo(xyz)
#'
#' @name sphere.convert
#' @concept sphere
#' @rdname sphere.convert
NULL

#' @rdname sphere.convert
#' @export
sphere.geo2xyz <- function(lat, lon){
xlat = as.double(lat)*pi/180
xlon = as.double(lon)*pi/180

x = cos(xlat)*cos(xlon)
y = cos(xlat)*sin(xlon)
z = sin(xlat)

outvec = c(x,y,z)
return(outvec/sqrt(sum(outvec^2)))
}

#' @rdname sphere.convert
#' @export
sphere.xyz2geo <- function(xyz){
x = as.double(xyz[1])
y = as.double(xyz[2])
z = as.double(xyz[3])

latitude  = 180*asin(z)/pi
longitude = 180*atan2(y,x)/pi

output = c(latitude, longitude)
return(output)
}


## Try the Riemann package in your browser

Any scripts or data that you put into this service are public.

Riemann documentation built on March 18, 2022, 7:55 p.m.