Computing the Oja Median in R : The Package OjaNP
Full text
JSS Journal of Statistical Software February 2020, Volume 92, Issue 8. doi: 10.18637/jss.v092.i08 Computing the Oja Median in R: The Package OjaNP Daniel Fischer Natural Resources Institute Finland & University of Tampere Karl Mosler University of Cologne Jyrki Möttönen University of Helsinki Klaus Nordhausen Vienna University of Technology Oleksii Pokotylo University of Cologne Daniel Vogel University of Aberdeen Abstract The Oja median is one of several extensions of the univariate median to the multivariate case. It has many desirable properties, but is computationally demanding. In this paper, we first review the properties of the Oja median and compare it to other multivariate medians. Then, we discuss four algorithms to compute the Oja median, which are implemented in our Rpackage OjaNP. Besides these algorithms, the package contains also functions to compute Oja signs, Oja signed ranks, Oja ranks, and the related scatter concepts. To illustrate their use, the corresponding multivariate oneand C-sample location tests are implemented. Keywords: Oja median, Oja signs, Oja signed ranks, Oja ranks, R,C++. 1. Introduction The univariate median is a popular location estimator. It is, however, not straightforward to generalize it to the multivariate case since no generalization is known that retains all properties of the univariate estimator, and therefore different generalizations emphasize different properties of the univariate median. So besides the Oja median described here, there are several other multivariate median concepts. Hayford (1902) suggested the first generalization by simply using the vector of the marginal medians. Other popular multivariate medians are Tukey’s median (Tukey 1975) and the spatial median (also known as L1median). The spatial median was initially defined as a bivariate median (Weber 1909,1929) and subsequently
2OjaNP: Computing the Oja Median in R extended to the general multivariate case. These and more multidimensional medians are surveyed in Small (1990) and Oja (2013). While the vector of marginal medians is quite easy to compute, the other multivariate medians are more computationally expensive. Particularly the Oja median (Oja 1983) has, despite its compelling statistical properties, not been used very often in practice so far, since it is difficult to compute. The main topic of this paper is to describe the R(RCore Team 2019) package OjaNP (Fischer, Mosler, Möttönen, Nordhausen, Pokotylo, and Vogel 2020), which provides several algorithms for the computation of the Oja median in Rand is available from the Comprehensive RArchive Network (CRAN) at https://CRAN.R-project.org/package=OjaNP. The outline of this paper is as follows. In Section 2.1, we review and compare some multivariate medians and show which properties of the univariate median is generalized by which multivariate median. Our main focus is on the Oja median, whose basic properties are discussed in Section 2.2, followed by an introduction to Oja signs and ranks (Section 2.3), Oja signed ranks (Section 2.4) and Oja sign and rank covariance matrices (Section 2.5). To demonstrate the application of the Oja median and its sign and rank concepts, Section 2.6 discusses oneand C-sample tests of location. In Section 3, we focus on the different algorithms provided by the package OjaNP to calculate the Oja median. Four different algorithms are available: two exact algorithms and two approximate algorithms based on different designs. Section 4shows how to use the package OjaNP in order to calculate the Oja median and related statistics. We provide simple examples, and additional benchmarks are calculated to analyze the performance of the implementations. A concept frequently encountered in this paper is affine equivariance. We use it in the sense of full-rank affine equivariance, which is common in robust statistics. Given the k-dimensional sample x1,...,xnwe let X= (x1. . . xn)>be the data matrix of dimension n×k, containing the data points as rows. Subsequently the data sample is identified with X. For a given affinelinear transformation T:Rk→Rk,x7→ Ax +bwith b∈Rkand A∈Rk×knon-singular, the data matrix Yof the transformed data T(x1), . . . , T(xn)is given by Y=T(X) = XA>+1b>, where 1denotes the n×1vector consisting of ones. We call an Rk-valued location statistic µ(X)affine equivariant, if µ(T(X)) = T(µ(X)) for all Tas above. (1) This applies analogously to set-valued location statistics µ, such as median sets. For a matrix-valued scatter statistic Staking on values in Rk×k, affine equivariance is commonly understood as S(T(X)) = AS(X)A>. For a more detailed introduction to affine equivariance, see Oja (2010).
Journal of Statistical Software 3 2. Oja median and related concepts 2.1. Oja median and other multivariate medians We start by introducing the univariate median for distributions. Given a distribution function F, let F−11 2−= inf x∈R:F(x)≥1 2and F−11 2+= sup x∈R:F(x)≤1 2. Then the median (set) of Fis given by the interval Med(F) = F−11 2−, F−11 2+.(2) Any point of the interval divides the distribution in two halves of equal probability weight and can represent the median. In case a unique selection is needed, we use the gravity center of the median set as a (single-point) median and denote it by the lower case symbol, med(F) = F−11 2++F−11 2− 2.(3) For a given sample X= (x1, . . . , xn), the median med(X)is obtained as med(X) = x(n+1 2)if nis odd , 1 2x(n 2)+x(n 2+1)if nis even , with x(i)being the ith order statistic. The latter definition is a special case of (3), obtained by taking Fto be the empirical distribution of X, which gives equal probability mass to each of the points x1, . . . , xn. Note that med(X)is an affine equivariant location statistic. Now let X= (X1,...,Xk)be a k-dimensional data set, where Xidenotes the ith column of Xand corresponds hence to the ith variable. The many existing notions of a k-variate median for such data have in common that they reduce to the univariate median for k= 1. Multivariate medians are generally non-unique, and we select, as above, the gravity center of the median set to obtain a unique representation. To the best of our knowledge the first generalization of the univariate median to the multivariate case is the vector of marginal medians mmed described in Hayford (1902). Definition 1. The vector of marginal medians mmed of the sample Xis defined as mmed(X) = (med(X1), . . . , med(Xp))>. The vector of marginal medians is easily computed but not affine equivariant. Example 1. A simple rotation of two-dimensional data visualizes the problem of not affine equivariant transformations. In the left part of Figure 1a two-dimensional data set is plotted and rotated around the center point +. At the same time the marginal median mmed is
4OjaNP: Computing the Oja Median in R ● ● ● ● ● ● ● ● ● í í í í í í í í ; < ● ● ● ● ● ● ● ; < ● ● ● ● ● ● ● ● í í Figure 1: Transformation of different multivariate medians. continuously plotted in red. An affine equivariant median would follow the rotation on a circle, meaning that it would not be affected by the degree of the rotation of the original data. The blue circle represents the Oja median and hence, illustrates the behavior of an affine equivariant median. In the right part of Figure 1we visualize the behavior of a rotation invariant estimator that is not scale invariant. We use the spatial median smed as an example for a not scale invariant statistic, see Definition 3below. For five example points the median is calculated and indicated as the center point, connected with lines to the data points. Then, the scale of those five points was transformed into the star shape on the right side of this plot. The location of the spatial median of the transformed data is indicated at the top of curved arrow. An affine equivariant median, however, would still be equivalent to the transformed median of the original data. The Oja median (and any other affine equivariant statistic) fulfills this criterion. The second generalization of the univariate median reviewed here is based on the fact that for a given sample Xand its median med(x)the equation n i=1 11(−∞,med(x)](xi)= n i=1 11[med(x),∞)(xi) holds, where 11A(xi)=1if xi∈Aand =0otherwise. This means, there are as many observations smaller than med(x)as are bigger, as it was pointed out, e.g., in Hotelling (1929). Given a k-dimensional data set, we consider halfspaces in every direction v∈Sk−1= {x|x=1}.LetHvdenote the “minimal” halfspace with normal vector vthat contains at least half of the data points, i.e., for any other halfspace ˜ Hvwith these properties Hv⊂˜ Hv holds. The intersection of all Hv,v∈Sk−1, forms the Tukey median Tmed(X). If the data are in general position, each such halfspace is bordered by a hyperplane through exactly k data points, and the Tukey median is in general no singleton. (A set of k-variate data is in general position if at most kof them lie on the same hyperplane.) The unique Tukey median tmed(X) is defined as the gravity point of this median set; see Tukey (1975)andDonoho and Gasko (1992). It is affine equivariant (Chen 1995) and can be introduced as the maximizer of a depth function as follows. Definition 2. Let Xbe a k-dimensional sample as above. For any x∈Rk, tdepX(x)= 1 nmin v=1 #{i:vxi≥vx}(4)
Journal of Statistical Software 5 is called the Tukey depth or location depth of xw.r.t. X. (Here #{S}denotes the cardinality of a set S.) The Tukey median Tmed(X)of the sample Xis defined as the maximizer of the Tukey depth, Tmed(X) = argmax x∈Rk{tdepX(x)}.(5) Another way to generalize the univariate median is to transfer its minimizing feature into higher dimensions. Consider again a univariate sample X= (x1, . . . , xn)in R. When minimizing the sum Pn i=1 |xi−x|over x∈Rwe obtain med(X) = argmin x∈R n X i=1 |xi−x|.(6) This minimizing feature of the univariate median is interpretable in two ways: It is the sum of absolute deviations, but it can also be viewed as the sum of one-dimensional simplices. The first interpretation will lead us to the spatial median, the second to the Oja median. The spatial median is presumably the most popular multivariate median and almost as old as the marginal median. Weber (1909) first described the spatial median and used it to solve an economic problem: He was looking for the best place of a distribution center in the sense, that the sum of distances between outposts and the distribution center becomes minimal. The k-dimensional extension of this approach is the spatial median. Definition 3. The spatial median smed of a k-variate data sample Xis defined as smed(X) = argmin x∈Rk(n X i=1 kxi−xk2).(7) Here k·k2denotes the Euclidean norm. The spatial median is sometimes also referred to as the L1median since it minimizes the L1norm of the n-variate vector of distances kxi−xk2. The spatial median is not affine equivariant as is visualized in the right part of Figure 1. Given an affine transformation T:R2→R2which transforms the data points from the left star-shaped figure into the right star-shaped figure, the center of each star is the spatial median of the data points. However, due to the lack of affine equivariance, the transformed spatial median T(smed(X)) (marked with the arrow) does not coincide with the spatial median smed(T(X)) of the transformed data set. The spatial median is rotation invariant, but not scale invariant. Oja (1983) introduced a multivariate median based on volumes of simplices. A k-dimensional simplex is the convex hull of (k+ 1) spanning points in general location. Let x1= (x1,1, . . . , x1,k)>,x2= (x2,1, . . . , x2,k)>,...,xk+1 = (xk+1,1, . . . , xk+1,k)> be k+ 1 points in general location from Rk. The volume V(x1,...,xk+1)of the simplex spanned by the points x1,...,xk+1 is then given by V(x1,...,xk+1) = abs 1 k!det 1 1 ··· 1 x1,1x2,1··· xk+1,1 x1,2x2,2··· xk+1,2 . . .. . .. . . x1,k x2,k ··· xk+1,k ,(8)
6OjaNP: Computing the Oja Median in R see, e.g., Stein (1966). Let Fbe a distribution on Rkhaving finite first moment. The Oja median Omed(F)of Fis defined as follows: for i.i.d. random vectors X1,...,Xk, each distributed with F, define the Oja median set as Omed(F) = argmin µ∈Rk E(V(X1,...,Xk,µ)). For data we have the following definition. Let Pn,k ={p= (i1, . . . , ik)|1≤i1<··· < ik≤n}(9) be the set of all ordered k-tuples out of {1, . . . , n},1≤k≤n. Definition 4. The Oja median Omed of a k-dimensional sample X= (x1. . . xn)>of size n>kis defined as Omed(X) = argmin x∈Rk X (i1,...,ik)∈Pn,k V(xi1,...,xik,x) .(10) This means that the Oja median of a k-dimensional sample is any point x∈Rkfor which the sum of simplex volumes over all combinations of possible kdata points is minimal. Note that Omed(X)equals Omed(FX), where FXis the empirical distribution on x1. . . xn. A unique version of the Oja median, denoted omed(X), is obtained by selecting the point of gravity of Omed(X). 2.2. Properties of the Oja median As other multivariate extensions of the median, the Oja median is not unique. In the bivariate case, we have the following result. Theorem 1. If nis even and the data points x1,...,xn∈R2are in general position, then the Oja median is unique. For details see Niinimaa (1995). The author also identifies a necessary and sufficient condition for the bivariate Oja median to be unique if nis odd. There appears to be no result for higher dimensions, but Oja (1999) conjectures that the Oja median is unique for even (odd) sample sizes if the dimension kis even (odd). Figure 2visualizes the Oja median in several smallsample data situations. Theorem 2. The Oja median of a sample is a convex set. See Oja and Niinimaa (1985) for details. Theorem 2is illustrated in Figure 3. Contour lines of the objective function (10) are plotted for two small data clouds. The objective function is convex for, both, even and odd sample sizes. Theorem 3. The Oja median is affine equivariant. Hence the Oja median is a proper location statistic in the sense of (1), i.e., omed(T(X)) = T(omed(X)). The next result can be found, e.g., in Arcones, Chen, and Gine (1994), Shen (2008).
Journal of Statistical Software 7 • • • Q • Q • • • • Q • • • • • Q • • • • • • Q • • • • • • • Q • • • • • • • Figure 2: Example plots for the bivariate case. • • • • • • • • • • • Figure 3: Contour plot for the bivariate case. Theorem 4. Under mild regularity conditions (including the existence of first moments) for an i.i.d. sample X=(x1...xn)from the distribution F, we have √n(omed(X)−omed(F)) →dNk(0,A−1B(A−1)), where Bis the Oja sign covariance matrix (OSCM) defined later in Section 2.5 at Fand Ais the expected covariance matrix between the Oja signs (osgn, see Section 2.3) and the optimal location score of F.
8OjaNP: Computing the Oja Median in R Median Affine equivariance Breakdown point Marginal median No 1/2 Tukey median Yes 1/(k+ 1) Spatial median No 1/2 Oja median Yes 0 Influence function Asymptotic distribution Marginal medians Niinimaa and Oja (1995)Babu and Rao (1988) Tukey median Romanazzi (2001)Bai and He (1999) Spatial median Niinimaa and Oja (1995)Möttönen, Nordhausen, and Oja (2010) Oja median Niinimaa and Oja (1995)Shen (2008) Table 1: Main properties of different multivariate medians. Concerning robustness the Oja median has the following properties. Theorem 5. Let Fhave a finite first moment. Then the Oja median has a bounded influence function. The influence function is for example given in Niinimaa and Oja (1995). The Oja median has an asymptotic breakdown point of 0. Niinimaa, Oja, and Tableman (1990) show the following: Theorem 6. The finite-sample breakdown-point of the bivariate Oja median is 2/(n+ 2). Table 1summarizes the main properties of the Oja median and other common multivariate medians discussed here. For another recent discussion of the different medians see also Oja (2013). 2.3. Oja signs and ranks Closely related to the median is the concept of signs and ranks. In this section we introduce the multivariate Oja sign and Oja rank and relate them to other multivariate signs and ranks, corresponding to the marginal and the spatial median. Again, to motivate the subsequent derivations we take a brief look at the univariate case. Suppose we have a univariate data set X= (x1. . . xn)>,n∈N. We call sgnX(x) = sgn(x−med(X)), x ∈R,(11) the sign of xw.r.t. the data sample X, where sgn is the univariate sign function (sgn(x) = x |x| if x6= 0 and zero otherwise), and med(X)is the univariate median of the sample X. There are several possibilities of suitably assigning ranks to the data points. By rnkX(x) = 1 n n X i=1 sgn(x−xi), x ∈R,(12) we define normalized central ranks, which may take on 2n+ 1 possible values ranging from −1to 1. In the following we call rnkX(x)simply the rank of xw.r.t. X.
Journal of Statistical Software 9 The median appropriately centers the data, i.e., the signs of the data points, centered by the median, sum up to zero: n X i=1 sgnX(xi) = n X i=1 sgn(xi−med(X)) = 0.(13) The sample mean does the same for the data points themselves. In other words we may say that the median has central rank zero, rnkX(med(X)) = −1 n n X i=1 sgn(xi−med(X)) = 0 ,(14) and, in this respect, is the most central point. Identities (13) and (14) (which are different formulations of the fact that half of the data lies above and below the median) provide the essential link between signs and ranks and the median and are a motivating principle behind the multivariate sign and rank functions we will introduce next. Note that a k-variate sign function should be a vector that can point in any direction of the k-dimensional space. The same holds for a multivariate rank function based on signs. An obvious extension of (11) to the multivariate setting is its componentwise application, leading to the marginal sign function. We call msgnX(x) = msgn(x−mmed(X)) the marginal sign of x∈Rkw.r.t. the k-variate data sample X= (x1. . . xn)>, where x= (x1. . . xk)>, msgn(x) = (sgn(x1). . . sgn(xk))>, and mmed(X)is the marginal median of X. An equally straightforward generalization is the spatial sign of xw.r.t. X: ssgnX(x) = ssgn(x−smed(X)),x∈Rk, where ssgn(x) = 1 kxkxif x6=0, 0if x=0, and smed(X)is the spatial median of X. The corresponding rank functions, the marginal rank mrnkXand the spatial rank srnkX, are obtained by replacing sgn in (12) by msgn and ssgn, respectively. The Oja sign is defined as follows. For 0≤k≤n, let Nn,k =n kand Pn,k as in (9). We call osgnX(x;m) = 1 Nn,k−1X (i1,...,ik−1)∈Pn,k−1∇xdet(xi1−m,...,xik−1−m,x−m)(15) =1 Nn,k−1X (i1,...,ik−1)∈Pn,k−1∇x det 1 1 . . . 1 1 m xi1. . . xik−1x! the Oja sign of the point x∈Rkw.r.t. the data sample Xand the center location mand osgnX(x) = osgnX(x;omed(X)) (16)
16 OjaNP: Computing the Oja Median in R is there to view the signs, signed ranks and ranks as scores, which replace the observations in the classical multivariate procedures. In principle, robust counterparts of any multivariate method can be derived this way. We demonstrate here the multivariate one-sample and Csample location tests. For that purpose, denote for a sample point xithe corresponding score s(xi;m), where mis an optional location w.r.t. to which the score is computed. The one-sample tests Assume X= (x1,...,xn)is a sample of size nfrom a k-variate symmetric distribution Fwith symmetry center µ. We are interested in testing the null hypothesis H0:µ=µ0against H1:µ6=µ0. Denote ¯ s=1 nPn i=1 s(xi,µ0)as the average of the score values under the null hypothesis and Σs=1 nPn i=1 s(xi,µ0)s(xi,µ0)>. The test statistic is then Q=n¯ s>Σ−1 s¯ s. Using Oja signs or Oja signed ranks as scores, this yields a straightforward extension of Hotelling’s classical one-sample T2-test. The test is invariant under affine transformations and asymptotically distribution-free and has a limiting χ2 kdistribution. Test decisions can also be based on permutation principles by randomly changing the signs of the scores. The tests are described in detail in Hettmansperger, Nyblom, and Oja (1994) and Hettmansperger, Möttönen, and Oja (1997). Similar tests can also be naturally constructed using marginal or spatial signs and signed ranks. See Puri and Sen (1971) and Oja (2010) for details. The C-sample tests Let Xc= (xc,1,...,xc,nc),c= 1, . . . , C, correspond to k-variate samples coming from C≥2 groups having distributions Fcthat differ only in location parameters µc,c= 1, . . . , C. The null hypothesis is µ1=. . . =µC, i.e., the Cgroups have the same location. Denote X= (X1,...,XC)as the combined sample and n=PC i=1 ni. Then ¯ sc=1 ncPnc j=1 s(xc,j), c= 1, . . . , C, is the average score value of group ccomputed w.r.t. to the location of the combined sample. Similarly, Σs=1 nPn i=1 s(xi)s(xi)>is computed for the combined groups. The test statistic is obtained as Q= C X c=1 nc¯ s> cΣ−1 s¯ sc. When Oja signs and Oja ranks are used as scores, C-sample tests for multivariate location are obtained that are asymptotically distribution-free and affine invariant. The limiting distribution of the test statistic is χ2 k(C−1), but pvalues can be obtained by permuting observations between the groups. The two tests are described in Hettmansperger and Oja (1994); Hettmansperger, Möttönnen, and Oja (1998). Similar tests based on other concepts of signs and ranks are described also in Puri and Sen (1971) and Oja (2010). 3. Description of the algorithms The package OjaNP contains four different algorithms to calculate the Oja median. Two exact algorithms and two approximate algorithms. The first exact algorithm was developed
Journal of Statistical Software 17 in Ronkainen, Oja, and Orponen (2003) as well as one of the approximate algorithms. The second exact algorithm (Mosler and Pokotylo 2015) is based on the first one: It accelerates the computation considerably by introducing bounds to the region of search. The numerical calculation is a non-trivial problem which consumes enormous calculation resources and hence, the exact algorithms are limited to small data situations only and, as a consequence, approximate algorithms are needed. These offer parameters to regulate the speed vs. accuracy trade-off, and the user has to decide from case to case which algorithm to choose with which tuning parameters. In Section 4we will give an overview over the several options and their effect onto the calculation precision and time. Before that, we are going to describe the four algorithms. 3.1. Exact algorithm Ronkainen et al. (2003) implemented the ideas from Niinimaa, Oja, and Nyblom (1992) and generalized them into higher dimensions based on the result described in Hettmansperger, Möttönen, and Oja (1999), whereby the vertices of the Oja median set are always located on intersections of hyperplanes that are spanned by data points. Ronkainen et al. (2003) constructed a Las Vegas algorithm as follows. (This is a simplified version; a more detailed description can be found in the original paper.) 1. Let Hbe the set of all (k−1)-dimensional hyperplanes spanned by the points in X. 2. Take the data point xc∈Xclosest to the mean as an initial candidate point. 3. Sample k−1hyperplanes out of Hsuch that the candidate point is on their intersection L. 4. Calculate the Oja depth of each intersection point between Land the hyperplanes in H. 5. Take the point x∗ cwith the highest Oja depth as next candidate point for the Oja median. 6. Repeat steps 3 to 5 until no improvement in the objective function is possible (or latest after nrepetitions). 7. The result for the exact Oja median is the last candidate point x∗ c. Ronkainen et al. (2003) focused on computational stability rather than efficiency. In case a candidate point is a data point, there are k−1possible intersection hyperplanes L. Instead of only following the one determined by the gradient of the objective function, the algorithm tries all possible ones. This algorithm finds just one of the vertices of the median set. While searching for the median, the algorithm may pass through several vertices of the median set, although it is not guaranteed that it visits all of them. The reason is that on step 5 only the first of possibly two points having highest Oja depth is taken as x∗ c. However, in case of a non-unique median, there exist two such points lying on an edge of the median set. To deliver all vertices of the median set, the algorithm can be modified as follows: It has to store both points as vertices and, in addition, check all lines passing through them.
18 OjaNP: Computing the Oja Median in R The algorithm is implemented in C++ and was initially published as stand-alone software on the personal webpage of Tommi Ronkainen but cannot be accessed anymore in that form. The implementation was modified such that it can be used directly from R. 3.2. Exact bounded algorithm Based on the exact algorithm of Ronkainen et al. (2003), Mosler and Pokotylo (2015) developed a faster exact algorithm. This algorithm uses the centered rank functions to build bounded regions which contain the median. The negative rank function −ornkX(x)is a vector that points in a direction of ascent of the depth function. It defines a hyperplane through x, on the positive side of which the Oja median is found. The halfspaces defined by the negative rank function are used to build a bounded region that contains the median. In this algorithm, these halfspaces are selected in an iterative way and the further search is restricted to their intersection. The hyperplanes bordering such a search region will be called bounding hyperplanes or simply bounds. The steps of the algorithm are as follows. (A more detailed description can be found in the original paper): 1. Let Hbe the set of all hyperplanes spanned by the points in X. 2. Create the initial rectangular bounded region B, limited by hyperplanes that are perpendicular to the coordinate axes and go through the maximal and minimal coordinates of the data points on these axes. 3. Iteratively reduce the bounded region Bby adding hyperplanes that go through a properly chosen central point of the region and have their normal vectors equal to the corresponding negative rank function. Specifically, the mean value of the bounds’ intersection points is selected as a central point in our implementation. The bounds of Bare cut off by newly added hyperplanes. 4. The region is reduced until the desired final volume of the bounded region is reached. 5. Add the bounds from Bto H. •At the first iteration: Take a random initial line Lon the border of B. •At further iterations: Sample k−1hyperplanes out of Hin that way, that the candidate point is on their intersection line L. 6. Calculate the Oja depth on each intersection point between Land hyperplanes in H that lies in the bounded region. 7. Take the point x∗ cwith the highest Oja depth as next candidate point for the Oja median. 8. Repeat steps 6 to 8 until no improvement in the objective function is possible (or latest after nrepetitions). 9. The result for the exact Oja median is the last candidate point x∗ c.
Journal of Statistical Software 19 0 5 10 15 20 25 30 0 50 100 150 200 250 Remaining Volume Time T_Total T_bounds T_count Figure 7: Dependence of calculation time on the size of the bounded region; time needed for bounding (Tbounds), for minimizing (Tcount), and total time (Ttotal). For comparison: total time of the first exact algorithm is about 340 seconds. The bounded regions reduce the complexity of the searching procedure by reducing the number of hyperplanes that cross the searching lines as well as the number of their intersections to be considered in the minimization procedure. The algorithm is driven by the desired final volume of the bounded region, which is the volume of the minimal rectangle containing the region and having edges parallel to the coordinate axes. As this parameter is reduced, the time needed to build the bounded regions (Tbounds) increases, while the minimization time (Tcount) decreases along with the number of hyperplanes and their intersections; see Figure 7. The total computing time (Ttotal) decreases rapidly with the volume, but then slowly grows again. Beyond some point the procedure becomes less efficient. For comparison, the original algorithm needs a total time Ttotal of ca. 340 seconds in this example. It appears that the fastest computation is obtained if bounds are imposed until the volume of the bounded region ranges around 10−8of the original volume. Note that the bounds may cut off some of the vertices of the median set. Moreover, if the central point of the bounded region lies in the median set on step 3, its negative rank function is zero, and this point is directly returned as a median, as on Figure 9. Due to limitations in computing memory and long calculation time, the exact algorithm and its bounded version are only able to calculate the Oja median for small data sets in low dimensions. For example, the calculation of the median in a data set of size 100 ×5needs 12 GB RAM. Therefore, approximate algorithms are needed. Obviously, the bounded algorithm can be stopped at any iteration and some mean value of the last bounded region be taken as an approximation of the Oja median. However, unlike the approximate algorithms presented below, this approach requires the calculation of all hyperplanes. Due to its high computational requirements it is less suited as an approximate algorithm for big data sets in high dimensions. For calculation speed reasons, also this algorithm is implemented in C++.
20 OjaNP: Computing the Oja Median in R 3.3. Grid-based algorithm The third algorithm, which calculates an approximation to the Oja median, was also proposed in Ronkainen et al. (2003). Technically it is a Monte Carlo algorithm. The algorithm lays a uniform grid over the data set. At each grid point a test is performed whether the point is a possible candidate. The amount of candidate points is reduced as long as only one grid point is left. This point is afterwards the center for a smaller but denser grid, where again each grid point is tested. The algorithm stops when the distance between two grid points gets smaller than a predefined parameter. A second tuning parameter is the significance level of the point tests. The steps of the algorithm are: 1. Create a grid Gwith equidistant knot distance h, covering the whole data set. 2. Choose randomly a set of hyperplanes, build the test statistic and test each of the grid knots in Gwhether it is an Oja median. 3. Remove those grid knots which have been tested not to be an Oja median. 4. If there is more than one knot left, sample additional hyperplanes and repeat the test for the remaining knots. 5. Repeat these steps until only one grid knot is left over. If the last test removes all remaining ones, take the last set. 6. Build a new grid around the last remaining old grid knot with equidistant knot distance h/2. 7. Repeat all these steps until the grid distance reaches a predefined threshold. 8. The last point is taken as the Oja median. For further details, especially about the testing procedure, we refer to the original paper Ronkainen et al. (2003). It may happen that there is continually more than one point left on step 5. In this case sampling of the additional hyperplanes may not help and the algorithm hangs. We restrict the algorithm to 5000 iterations on step 5, after which the grid point with the best test statistic is passed to step 6 and the number of iterations is reduced to 100. We repeat until the grid threshold is reached and return the grid point with the best test statistic as an Oja median approximation. The grid-based algorithm is part of the initially stand-alone tool that contains also the implementation of the first exact algorithm. It is as well written in C++ to speed up the computational burden. 3.4. Evolutionary algorithm The fourth algorithm to calculate the Oja median is an evolutionary algorithm, which is based on mutations of the latest candidate points. It was developed by the Department of Computer Science, Efficient Algorithms and Complexity Theory at the TU Dortmund, but has not been published before. The algorithm works as follows: 1. Set the level of initial mutation variance σ2 0.
Journal of Statistical Software 21 2. Take 10 randomly chosen observations from X. 3. Evaluate the objective function for all these points and take the minimum as starting candidate point η. 4. Choose k+ 1 random numbers x1, . . . , xk, l from a N(0, σ2 0)distribution and calculate the k-variate mutation vector ν=|l| qx2 1+···+x2 k (x1, . . . , xk)>. The mutation vector νhas a normally distributed length with given variance σ2 0and uniformly distributed direction. 5. Calculate mmutation points ηi=η+νifor i= 1, . . . , m from the last candidate point η. 6. Calculate the ratio rhow often the objective function is bigger at the mutations than at η. 7. If r > 0.2then σ2 0·κelse σ2 0·1 κfor some κ > 1. 8. Choose as a new candidate point the mutation with the smallest objective function value. 9. Repeat steps 4 to 8 until the variance for the next mutation drops under a predefined value s. 10. If the algorithm has not terminated after ntsteps, stop the calculation. Step 7 controls the dynamic of the mutation. If at more than 20% of the mutations the objective function has a smaller value than at the last candidate point, the algorithm increases the variability of the mutation; hence the search area is enlarged. Step 10 ensures that the algorithm terminates in any case. As the other implementations also this algorithm is implemented in C++. 3.5. Other algorithms for the Oja median in R There are other implementations of algorithms for the Oja median available in R, but they are mostly restricted to two dimensions. The function med in the package depth (Genest, Masse, and Plante 2019) uses the Fortran code of Niinimaa et al. (1992) and is restricted to the bivariate case. Another method to compute the Oja median was suggested by Roger Koenker on R-help on 2003-08-16 (see https://hypatia.math.ethz.ch/pipermail/r-help/2003-August/037702. html using the quantreg package; Koenker 2019): R> oja.median <- function(x) { + n <- dim(x)[1] + A <- matrix(rep(1:n, n), n)
22 OjaNP: Computing the Oja Median in R + i <- A[col(A) < row(A)] + j <- A[n + 1. - col(A) > row(A)] + xx <- cbind(x[i, ], x[j, ]) + y <- xx[, 1] * xx[, 4] - xx[, 2] * xx[, 3] + z1 <- (xx[, 4] - xx[, 2]) + z2 <- -(xx[, 3] - xx[, 1]) + return(quantreg::rq(y ~ cbind(z1, z2) - 1)$coef) + } 4. The Rpackage OjaNP The main purpose of the OjaNP package is to provide users with the possibility to compute the Oja median. The package includes, however, also other useful functions. The main functions of the package are visualized in Figure 8. Most of the function names are self-explanatory. For details about them we refer to the corresponding help pages. In the following we will first explain the function ojaMedian and its options in detail. Then we demonstrate the use of some of the functions with a small but illustrative data set, which is also contained in the package. ojaMedian ojaSign ojaSignRank ojaRank ojaSCM ojaRCM oja1sampleTest ojaCsampleTest Figure 8: The main functions in package OjaNP.
Journal of Statistical Software 23 • • • • • • • ; < ••••• • • • • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••••• • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••• • • • • • • • • • • • ••• • • • •• •••• • • • • • • • • • • • • • •• • •• • • • • • •• • •• • • • •••• •• • •• • • • •• • •• • • • • •• • •• •• •• •••• •• • • • ••• • • • • • • • •• • •• • • • •• •• • • ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• • ; < • • • • • • • • •• • • •• • • •• • • •• •• •• • • •• •• ••• •• •• • •• •• •• •• • • •• • • • • • • • • • • • • • • •• • • ••••••• • •••• •• •• •• • •• • • • • • • • •• •• • • ••••• • • • • • • • • • • • • • •• • • • • •• • • • • • • • • • • • •• •• • • • • • • • • • • • • • • • • •• • • • • • • • • • • • •• • • ••••• • • • • • • • • • • • • •• • • • • • • • • ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• Figure 9: Example plot for all four algorithms: exact algorithms (red), evolutionary algorithm (blue), grid algorithm (green). The convex median set is marked with yellow. 4.1. The computation of the Oja median in OjaNP The main function of OjaNP is ojaMedian(X, alg = "evolutionary", sp = 1, na.action = na.fail, control = ojaMedianControl(...), ...) The user can choose via the alg option between four algorithms to calculate the Oja median. Furthermore, we have an option to calculate the Oja median repeatedly and average these results in order to receive less varying results. The amount of repetitions can be controlled with the sp parameter. In what follows we are going to explain the different parameters which control the flow of the different algorithms in detail and give insights how to choose the parameters in a given data situation. The default algorithm of the ojaMedian function is the evolutionary algorithm. The evolutionary algorithm relaxes the affine equivariance property of the Oja median. In order to restore it we first perform a scatter matrix transformation to obtain an invariant coordinate system (implemented in ICS;Nordhausen, Oja, and Tyler 2008), apply the algorithm to the transformed data and re-transform afterwards. That way we restore the affine equivariance for this implementation. Figure 9shows the exemplary outcome of the four implemented algorithms in simple data situations for the unique (left) and non-unique (right) case. We have chosen simple data situations with 6 and 7 data points, run the different algorithms 500 times and plotted the outcome into the figure. In the left part of Figure 9we have a data set that has a unique Oja median. This one is correctly determined by the exact implementations (red), whereas both approximate algorithms have a systematic behavior which does not differ strongly from the non-unique case in the right side of the graphic. The evolutionary algorithm (blue) determines the Oja median always along lines of intersection with the result of a bordered area. The right-hand side of Figure 9exhibits data that have a non-unique Oja median. Here the two exact algorithms find a vertex (red) of the median set, while the evolutionary algorithm (blue) yields any point of the border of the median set, that is the area with the lowest Oja depth (yellow). The grid algorithm, however, terminates in this case usually within the convex median set.
24 OjaNP: Computing the Oja Median in R • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ; < • • • • • •• •• • • • • ••• • • • •• • •• • • • • •• •• • • • •• • •• • •• •• • • •• • • • • • • • •• • •• •• • • •• • • • • • •• • •• • • • • •• • • • • • •• •• •• • •• • •• • ••• •• • • • •• • • • • • • •• •• ••• • • • •• • •• • • • • • • • •• • •• •• • •• • • •• • • • •• • ••• •• • • • • • • • • • • • • •• • • •• • • • • • •• • • •• • • • • •• •• • • •• • • •• • • •• • •• • • • •• •• • •• •• • •• • • • • •• • • • • • • • •• •• • • • •• • •• • • • • •• • •• • • • •• • • •• • •• •• • • • • • •• •• • • •• • • •• • • • •• •• • • • • • • • •• • • • •• • •• • •• •• • •• • •• • • • • • •• • • • •• •• • • •• • •• • • • • •• • • • • • •• • • •• • • • • • • • •• • •• • • • • • •• • • •• •• • • •• • • • • • • ••• • • • • •• •• •• • • •• •• •• • • •• •• • • •• • • •• • • • •• • • •• •• • •• •• • •• • •• • • • •• • • ••• •• • •• ••• • •• • ••• • • • • •• •• ••• • • • •••• • •• • • ••••••• • • • ••••••••• • ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• • • • ••••• • • • • • • • • • • • • • • • ••••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••••• • • • • • • • ••• • • • • • • • • • ••• •••••• • • • • • ••• • ••••••• • • • • • • • • • ••• • • •• •• • • • • • • • ••••• • ••••••••• • ••• • • • • • • • • • • • • • • • • • • • • • • • ••••••• • • • • • ••• • ••• • ••• • ••• • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • • ••••• • • • • • ••• • ••• • • • ••••••••••••• • • • ••• • ••••• • • • ••• •• •••• •• •• •••••••••••••••••••••••••••••••••••• • ••••••••••••••••••••••••••••••••• • •• •••••••• • •• • •• • • • • • •• • • • • • • •• •• • • •• • • • • • • • • • •• • • • • • • • • • • •• • • • •• • • • • •• •• •• • • •••• •• • • • •• •• • • • • • • • • • • •• • • • • • • • • • • • • • • • • ••• • • • • • • • •• • • •• •• • • • • • • • • • • • • • • • •• • • • • • • • • • • •• • • • •• • •• • •• • • • • • • •• • • • • •• •• • • • • • •• • • • • • • ••• •• • • • • • • • ; < • • • • • • •• • • • • • • • • • • • • • ••• • • •• • • • • • • • • • • • • •• • • • • • • • • • • •• • • • • •• • • • •• •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • •• • • • • • •• • •• • • • •• • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • •• • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • •• • • • • • • • • • • • ••• • • • • • • • • • • • • • • • •• • • •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • •• • • • • • • • • • • •• • • • • • • • • •• • • • •• • •• • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • •• • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • •• • • • •• • • •• • • • • • • • • • • • • ••• • •• • • • • • • • • • ••• • • • • • • • • • •• • • • • • • • • ••• ••••• • ••••• • ••••• • • • •••••••• • • • • • • • • • • ••••••••••••••••• • • • • • ••••• • ••• • • • • • • • • • ••• • •••••••••••• • • • ••••••••••• • • • • • • • ••••• • ••••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • ••• • ••• • • • • ••• • • ••• • • • • • • • • • • ••• • • • • • • • • • • • • • •••••• • • • • ••• • • • • • • • • • • • • ••• • ••• • • • • • ••• • • • • • • • • ••• • ••••••• • • • • • • • ••• •••••••• • •••••••• • • ••• • •• • • • • •• • • • • • • • • • • • ••••• • ••• • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 2MD([DFW (YR$OJ *ULG Figure 10: Example plot for all four algorithms: exact algorithms (red), evolutionary algorithm (blue), grid algorithm (green). In a next step we are going to analyze the outcome of the algorithms in more complex data situations. The first typical data situation is a multivariate normal distributed data cloud and we calculate the Oja median with the exact (red), the grid (green) and the evolutionary algorithm (blue). As we can see the evolutionary algorithm has the biggest variation within the results, but we have to keep in mind that the default setting of our function is tuned to calculate fast results. In the next paragraph we will discuss how to get more precise results in trade-off for calculation time. For most applications the faster, less accurate settings appear to be preferable. The second data set of interest consists of two data clouds, which are far away from each other. The motivation for that is to see whether the algorithms take this into account or if they get stuck in one cloud. As you can see in the right part of the figure, that all included algorithms are capable of calculating the Oja median to be in the center between both data clouds. Let us now compare the runtime of the two approximate algorithms depending on the dimension and size of the data in the bivariate case. We do not perform a runtime analysis for the exact algorithms. Their performance depends strongly on the computer used (particularly memory), much more so than the approximate ones. On an average computer, bivariate problems up to 1,000 observations are solvable in acceptable time, but for higher dimensions and sizes this decision has to be made from case to case. For example in the seven-dimensional case with tens of observations it takes minutes for the exact algorithm to find the solution. The biggest limitation of the exact algorithm is the memory required to store all hyperplanes. Even if we would allow infinite calculation time, the algorithm would still not be able to calculate the exact Oja median in more complex data situations because it cannot pre-calculate and access the total amount of hyperplanes, and hence we are facing a corresponding address space problem. This also applies to the exact bounded algorithm. Compared with the original exact algorithm, the bounded one finds the solution approximately two to five times faster. In order to analyze the runtime for different dimensions we simulated 10,000 multivariate normal distributed random numbers (with μ=0,Σ=I). The approximate grid algorithm (solid green line) is only able to calculate the Oja median up to 5 dimensions in acceptable time for this amount of data; this is why we did not take higher dimension situations into
Journal of Statistical Software 25 'LPHQVLRQVVL]HQ 7LPHLQ6HFRQGV 6L]HWLPHVA 7LPHLQ6HFRQGV • •••• • •• • • •• • • • •• • • • • • • • • • •• ••• •• •• • • •• • • • • •• • • • •• • •• • •• • • • • • ••• •• • • • • • • • • • •• • • • • • • • • • • • • • • •••• • •• •• • •••• •• •••• • • •• •• • •• • •• • • • •• • • • •• • • • • • • • • • • • • • • • •••••••• •• •••• • ••••• •• ••• ••• • ••• • • • • ••• • •• • •• • • • • • • • • • •• • • • •• • • • • • • • • • • • •• • • • • • • • • • ••• • • • • • • • • • • • • • • • • (YR$OJ (YR$OJZRFKHFNV *ULG Figure 11: Runtimes for the approximate algorithms: the approximate grid algorithm (solid green line), the evolutionary algorithm (solid blue line), the raw method without transformations or validation checks (dotted blue line). consideration. The evolutionary algorithm (solid blue) can calculate in the same time the Oja median of a 35-dimensional data set, and even higher dimensional problems are solvable. Since the ICS step (and there especially the data validation checks) takes a lot of computational time we implemented also a raw command to access the algorithms directly without transformations or validation checks (dotted blue line). In the analysis of different dimensions the raw algorithm does not bring huge advantages, but in the runtime analysis concerning the size of the data in the bivariate case (right-hand side of Figure 11)wecandetectahuge advantage for more than 105observation. The raw algorithm is even able to calculate the Oja median for sample size 5·107. Hence, our advice in time-critical situations with high sample size is to use the raw method. If the affine equivariant property is still required, we advise to perform beforehand ICS separately without performing the included data validation steps. The evolutionary algorithm has many tuning parameters, some of which control its accuracy. As we have seen in Figure 10, the default settings for these tuning parameters are preset to deliver fast results. In trade-off for higher computational time, the user can adjust the settings to obtain a more precise algorithm. The key parameters are useAllSubsets, nSubsetsUsed,andsigmaLog10Dec. The latter is the main abort criterion of the algorithm. It forces the algorithm to stop if the logarithmized initial variance differs more then the value of sigmaLog10Dec from the actual logarithmized variance. In other words, when the variance of the mutation vector is getting small enough, the algorithm stops. The settings for useAllSubsets and nSubsetsUsed control how many spanned hyperplanes should be taken into account during the calculation of the Oja median. Since the total amount of possible hyperplanes could be huge (it is n kfor nobservations in the k-variate case), the flag for useAllSubsets should be used carefully. It is more advisable to control this with the argument nSubsetsUsed. Raise this value together with sigmaLog10Dec for more precise values, lower them for faster results. The dynamics of the evolutionary algorithm are controlled via the parameters sigmaInit, sigmaAda and adaFactor. All these take control over the variance adjustments of the mutation vector. The parameter sigmaInit sets the initial variance of the mutation, the settings for sigmaAda control after how many mutation steps the mutation variance is adjusted and
32 OjaNP: Computing the Oja Median in R 5. Conclusions There are many different multivariate medians. In this paper we explained how the different medians extend different properties of the univariate median to the multivariate case. The Oja median has very convincing statistical properties, but is also among the computationally more challenging ones. We described the Rpackage OjaNP, which provides four different algorithms for its computation. Along with the concept of the Oja median comes the notion of Oja signs and ranks and multivariate scatter estimators based upon them. The package provides also functions for these useful tools, which can then be used for robust multivariate inferential procedures. As examples, we described and implemented the one-sample and the C-sample location test based on Oja signs and ranks and showed their practical use in simple and complex data situations. Acknowledgments The work of Klaus Nordhausen was supported by the Academy of Finland (grant 268703). Oleksii Pokotylo is supported by the Cologne Graduate School of Management, Economics and Social Sciences. The work of Daniel Vogel was supported by the DFG collaborate research grant SFB 823. The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources. References Arcones MA, Chen Z, Gine E (1994). “Estimators Related to U-Processes with Applications to Multivariate Medians: Asymptotic Normality.” The Annals of Statistics,22(3), 1460– 1477. doi:10.1214/aos/1176325637. Babu GJ, Rao CR (1988). “Joint Asymptotic Distribution of Marginal Quantiles and Quantile Functions in Samples from a Multivariate Population.” Journal of Multivariate Analysis, 27(1), 15–23. doi:10.1016/0047-259x(88)90112-1. Bai ZD, He X (1999). “Asymptotic Distributions of the Maximal Depth Estimators for Regression and Multivariate Location.” The Annals of Statistics,27(5), 1616–1637. doi: 10.1214/aos/1017939144. Chen Z (1995). “Robustness of the Half-Space Median.” Journal of Statistical Planning and Inference,46(2), 175–181. doi:10.1016/0378-3758(94)00105-5. Donoho D, Gasko M (1992). “Breakdown Properties of Location Estimates Based on HalfSpace Depth and Projected Outlyingness.” The Annals of Statistics,20(4), 1803–1827. doi:10.1214/aos/1176348890. Filzmoser P, Fritz H, Kalcher K (2018). pcaPP: Robust PCA by Projection Pursuit.Rpackage version 1.9-73, URL https://CRAN.R-project.org/package=pcaPP. Fischer D, Mosler K, Möttönen J, Nordhausen K, Pokotylo O, Vogel D (2020). OjaNP: Multivariate Methods Based on the Oja Median and Related Concepts.Rpackage version 1.0-0, URL https://CRAN.R-project.org/package=OjaNP.
Journal of Statistical Software 33 Genest M, Masse JC, Plante JF (2019). depth: Depth Functions Tools for Multivariate Analysis.Rpackage version 2.1-1.1, URL https://CRAN.R-project.org/package=depth. Hayford J (1902). “What Is the Center of an Area or the Center of a Population.” Journal of the American Statistical Association,8(58), 47–58. doi:10.2307/2276137. Hettmansperger TP, McKean JW (2011). Robust Nonparametric Statistical Methods. 2nd edition. CRC Press, Boca Raton. Hettmansperger TP, Möttönen J, Oja H (1997). “Affine-Invariant Multivariate One-Sample Signed-Rank Tests.” Journal of the American Statistical Association,92(440), 1591–1600. doi:10.1080/01621459.1997.10473681. Hettmansperger TP, Möttönen J, Oja H (1999). “The Geometry of the Affine Invariant Multivariate Sign and Rank Methods.” Journal of Nonparametric Statistics,11(1–3), 271– 285. doi:10.1080/10485259908832784. Hettmansperger TP, Möttönnen J, Oja H (1998). “Affine Invariant Multivariate Rank Tests for Several Samples.” Statistica Sinica,8, 785–800. Hettmansperger TP, Nyblom J, Oja H (1994). “Affine Invariant Multivariate One-Sample Sign Tests.” Journal of the Royal Statistical Society B,56(1), 221–234. doi:10.1111/j. 2517-6161.1994.tb01973.x. Hettmansperger TP, Oja H (1994). “Affine Invariant Multivariate Multisample Sign Tests.” Journal of the Royal Statistical Society B,56(1), 235–249. doi:10.1111/j.2517-6161. 1994.tb01974.x. Hotelling H (1929). “Stability in Competition.” The Economic Journal,39(153), 41–57. doi:10.2307/2224214. Koenker R (2019). quantreg: Quantile Regression.Rpackage version 5.54, URL http: //CRAN.R-project.org/package=quantreg. Mosler K, Pokotylo O (2015). “Computation of the Oja Median by Bounded Search.” In K Nordhausen, S Taskinen (eds.), Modern Nonparametric, Robust and Multivariate Methods, pp. 185–203. Springer-Verlag. Möttönen J, Nordhausen K, Oja H (2010). “Asymptotic Theory of the Spatial Median.” In J Antoch, M H˘usková, PK Sen (eds.), Nonparametrics and Robustness in Modern Statistical Inference and Time Series Analysis: A Festschrift in Honor of Professor Jana Jurecková, volume 7, pp. 182–193. Niinimaa A (1995). “Bivariate Generalizations of the Median.” In EM Tiit, T Kollo, H Niemi (eds.), Multivariate Statistics and Matrices in Statistics, pp. 163–180. VSP BV, Zeist. Niinimaa A, Oja H (1995). “On the Influence Functions of Certain Bivariate Medians.” Journal of the Royal Statistical Society B,57(3), 565–574. doi:10.1111/j.2517-6161. 1995.tb02048.x. Niinimaa A, Oja H, Nyblom J (1992). “Algorithm AS 277: The Oja Bivariate Median.” Journal of the Royal Statistical Society C,41(3), 611–633. doi:10.2307/2348099.
34 OjaNP: Computing the Oja Median in R Niinimaa A, Oja H, Tableman M (1990). “The Finite-Sample Breakdown Point of the Oja Bivariate Median and of the Corresponding Half-Samples Version.” Statistics & Probability Letters,10(4), 325–328. doi:10.1016/0167-7152(90)90050-h. Nordhausen K, Möttönen J, Oja H (2018a). MNM: Multivariate Nonparametric Methods. An Approach Based on Spatial Signs and Ranks.Rpackage version 1.0-3, URL https: //CRAN.R-project.org/package=MNM. Nordhausen K, Oja H (2011). “Multivariate L1Methods: The Package MNM.” Journal of Statistical Software,43(5), 1–28. doi:10.18637/jss.v043.i05. Nordhausen K, Oja H (2018). “Robust Nonparametric Inference.” Annual Review of Statistics and Its Application,5(1), 473–500. doi:10.1146/annurev-statistics-031017-100247. Nordhausen K, Oja H, Tyler DE (2008). “Tools for Exploring Multivariate Data: The Package ICS.” Journal of Statistical Software,28(6), 1–31. doi:10.18637/jss.v028.i06. Nordhausen K, Sirkiä S, Oja H, Tyler DE (2018b). ICSNP: Tools for Multivariate Nonparametrics.Rpackage version 1.1-1, URL https://CRAN.R-project.org/package=ICSNP. Oja H (1983). “Descriptive Statistics for Multivariate Distributions.” Statistics & Probability Letters,1(6), 327–332. doi:10.1016/0167-7152(83)90054-8. Oja H (1999). “Affine Invariant Multivariate Sign and Rank Tests.” Scandinavian Journal of Statistics,26(3), 319–343. doi:10.1111/1467-9469.00152. Oja H (2010). Multivariate Nonparametric Methods with R. An Approach Based on Spatial Signs and Ranks. Springer-Verlag, New York. Oja H (2013). “Multivariate Median.” In C Becker, R Fried, S Kuhnt (eds.), Robustness and Complex Data Structures. Festschrift in Honour of Ursula Gather, pp. 3–16. SpringerVerlag, Berlin. Oja H, Niinimaa A (1985). “Asymptotic Properties of the Generalized Median in the Case of Multivariate Normality.” Journal of the Royal Statistical Society B,47(2), 372–377. doi:10.1111/j.2517-6161.1985.tb01366.x. Ollila E, Croux C, Oja H (2004). “Influence Function and Asymptotic Efficiency of the Affine Equivariant Rank Covariance Matrix.” Statistica Sinica,14(1), 297–316. Ollila E, Oja H, Croux C (2003). “The Affine Equivariant Sign Covariance Matrix: Asymptotic Behavior and Efficiencies.” Journal of Multivariate Analysis,87(2), 328–355. doi:10.1016/ s0047-259x(03)00045-9. Puri ML, Sen PK (1971). Nonparametric Methods in Multivariate Analysis. John Wiley & Sons, New York. RCore Team (2019). R: A Language and Environment for Statistical Computing.RFoundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Romanazzi M (2001). “Influence Function of Halfspace Depth.” Journal of Multivariate Analysis,77(1), 138–161. doi:10.1006/jmva.2000.1929.
Journal of Statistical Software 35 Ronkainen T, Oja H, Orponen P (2003). “Computation of the Multivariate Oja Median.” In R Dutter, P Filzmoser, U Gather, PJ Rousseeuw (eds.), Developments in Robust Statistics: Proceedings of the International Conference on Robust Statistics (ICORS’01, Stift Vorau, Austria, July 2001), pp. 344–359. Springer-Verlag, Berlin Heidelberg. Shen G (2008). “Asymptotics of Oja Median Estimate.” Statistics & Probability Letters, 78(14), 2137–2141. doi:10.1016/j.spl.2008.02.004. Sievert C, Parmer C, Hocking T, Chamberlain S, Ram K, Corvellec M, Despouy P (2019). plotly: Create Interactive Web Graphics via plotly.js.Rpackage version 4.9.1, URL https: //CRAN.R-project.org/package=plotly. Small CG (1990). “A Survey of Multidimensional Medians.” International Statistical Review, 58(3), 263–277. doi:10.2307/1403809. Stein P (1966). “A Note on the Volume of a Simplex.” The American Mathematical Monthly, 73(3), 299–301. doi:10.2307/2315353. Tukey JW (1975). “Mathematics and the Picturing of Data.” In Proceedings of the International Congress of Mathematicians, volume 2, pp. 523–531. Vancouver. Visuri S, Koivunen V, Oja H (2000). “Sign and Rank Covariance Matrices.” Journal of Statistical Planning and Inference,91(2), 557–575. doi:10.1016/s0378-3758(00)00199-3. Visuri S, Ollila E, Koivunen V, Möttönen J, Oja H (2003). “Affine Equivariant Multivariate Rank Methods.” Journal of Statistical Planning and Inference,114(1–2), 161–185. doi: 10.1016/s0378-3758(02)00469-x. Vogel D, Fried R (2008). “Estimating Partial Correlations Using the Oja Sign Covariance Matrix.” In P Brito (ed.), COMPSTAT 2008 – Proceedings in Computational Statistics, volume II, pp. 721–729. Physica-Verlag, Heidelberg. Vogel D, Fried R (2011). “Elliptical Graphical Modelling.” Biometrika,98(4), 935–951. doi:10.1093/biomet/asr037. Vogel D, Köllmann C, Fried R (2008). “Partial Correlation Estimates Based on Signs.” In J Heikkonen (ed.), Proceedings of the 1st Workshop on Information Theoretic Methods in Science and Engineering. TICSP Series # 43. Weber A (1909). Über den Standort der Industrien. Mohr, Tübingen. Weber A (1929). Theory of the Location of Industries. The University of Chicago Press, Chicago.
36 OjaNP: Computing the Oja Median in R Affiliation: Daniel Fischer Natural Resources Institute Finland (Luke) Production Systems Jokioinen, Finland and School of Health Sciences University of Tampere Tampere, Finland E-mail: [email protected] Journal of Statistical Software http://www.jstatsoft.org/ published by the Foundation for Open Access Statistics http://www.foastat.org/ February 2020, Volume 92, Issue 8 Submitted: 2016-05-06 doi:10.18637/jss.v092.i08 Accepted: 2018-09-06