R version 3.6.0 Patched (2019-05-14 r76502) -- "Planting of a Tree" Copyright (C) 2019 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(cluster) > > tools::assertWarning(eh <- ellipsoidhull(cbind(x=1:4, y = 1:4)), verbose=TRUE) #singular Error in Fortran routine computing the spanning ellipsoid, probably collinear data Asserted warning: algorithm possibly not converged in 5000 iterations > eh ## center ok, shape "0 volume" --> Warning 'ellipsoid' in 2 dimensions: center = ( 2.5 2.5 ); squared ave.radius d^2 = 0 and shape matrix = x y x 1.25 1.25 y 1.25 1.25 hence, area = 0 ** Warning: ** the algorithm did not terminate reliably! most probably because of collinear data > stopifnot(volume(eh) == 0) > > set.seed(157) > for(n in 4:10) { ## n=2 and 3 still differ -- platform dependently! + cat("n = ",n,"\n") + x2 <- rnorm(n) + print(ellipsoidhull(cbind(1:n, x2))) + print(ellipsoidhull(cbind(1:n, x2, 4*x2 + rnorm(n)))) + } n = 4 'ellipsoid' in 2 dimensions: center = ( 2.66215 0.82086 ); squared ave.radius d^2 = 2 and shape matrix = x2 1.55901 0.91804 x2 0.91804 0.67732 hence, area = 2.9008 'ellipsoid' in 3 dimensions: center = ( 2.50000 0.74629 2.95583 ); squared ave.radius d^2 = 3 and shape matrix = x2 1.25000 0.72591 1.8427 x2 0.72591 0.52562 1.5159 1.84268 1.51588 4.7918 hence, volume = 3.1969 n = 5 'ellipsoid' in 2 dimensions: center = ( 3.0726 1.2307 ); squared ave.radius d^2 = 2 and shape matrix = x2 2.21414 0.45527 x2 0.45527 2.39853 hence, area = 14.194 'ellipsoid' in 3 dimensions: center = ( 2.7989 1.1654 4.6782 ); squared ave.radius d^2 = 3 and shape matrix = x2 1.92664 0.40109 1.4317 x2 0.40109 1.76625 6.9793 1.43170 6.97928 28.0530 hence, volume = 26.631 n = 6 'ellipsoid' in 2 dimensions: center = ( 3.04367 0.97016 ); squared ave.radius d^2 = 2 and shape matrix = x2 4.39182 0.30833 x2 0.30833 0.59967 hence, area = 10.011 'ellipsoid' in 3 dimensions: center = ( 3.3190 0.7678 3.2037 ); squared ave.radius d^2 = 3 and shape matrix = x2 2.786928 -0.044373 -1.1467 x2 -0.044373 0.559495 1.5496 -1.146728 1.549620 5.5025 hence, volume = 24.804 n = 7 'ellipsoid' in 2 dimensions: center = ( 3.98294 -0.16567 ); squared ave.radius d^2 = 2 and shape matrix = x2 4.62064 -0.83135 x2 -0.83135 0.37030 hence, area = 6.3453 'ellipsoid' in 3 dimensions: center = ( 4.24890 -0.25918 -0.76499 ); squared ave.radius d^2 = 3 and shape matrix = x2 4.6494 -0.93240 -4.0758 x2 -0.9324 0.39866 1.9725 -4.0758 1.97253 10.4366 hence, volume = 16.152 n = 8 'ellipsoid' in 2 dimensions: center = ( 3.6699 -0.4532 ); squared ave.radius d^2 = 2 and shape matrix = x2 9.4327 -2.5269 x2 -2.5269 3.7270 hence, area = 33.702 'ellipsoid' in 3 dimensions: center = ( 4.22030 -0.37953 -1.53922 ); squared ave.radius d^2 = 3 and shape matrix = x2 7.5211 -1.4804 -6.6587 x2 -1.4804 2.6972 11.8198 -6.6587 11.8198 52.6243 hence, volume = 84.024 n = 9 'ellipsoid' in 2 dimensions: center = ( 5.324396 -0.037779 ); squared ave.radius d^2 = 2 and shape matrix = x2 10.1098 -1.3708 x2 -1.3708 2.1341 hence, area = 27.885 'ellipsoid' in 3 dimensions: center = ( 5.44700 -0.12504 -1.13538 ); squared ave.radius d^2 = 3 and shape matrix = x2 7.0364 -1.2424 -5.5741 x2 -1.2424 1.7652 7.3654 -5.5741 7.3654 31.5558 hence, volume = 64.161 n = 10 'ellipsoid' in 2 dimensions: center = ( 4.85439 0.28401 ); squared ave.radius d^2 = 2 and shape matrix = x2 13.932 0.64900 x2 0.649 0.95132 hence, area = 22.508 'ellipsoid' in 3 dimensions: center = ( 5.12537 0.25024 0.86441 ); squared ave.radius d^2 = 3 and shape matrix = x2 9.29343 0.56973 1.4143 x2 0.56973 0.76519 1.8941 1.41427 1.89409 6.3803 hence, volume = 73.753 > > set.seed(1) > x <- rt(100, df = 4) > y <- 100 + 5 * x + rnorm(100) > ellipsoidhull(cbind(x,y)) 'ellipsoid' in 2 dimensions: center = ( -1.3874 93.0589 ); squared ave.radius d^2 = 2 and shape matrix = x y x 32.924 160.54 y 160.543 785.88 hence, area = 62.993 > z <- 10 - 8 * x + y + rnorm(100) > (e3 <- ellipsoidhull(cbind(x,y,z))) 'ellipsoid' in 3 dimensions: center = ( -0.71678 96.09950 111.61029 ); squared ave.radius d^2 = 3 and shape matrix = x y z x 26.005 126.41 -80.284 y 126.410 616.94 -387.459 z -80.284 -387.46 254.006 hence, volume = 301.25 > d3o <- cbind(x,y + rt(100,3), 2 * x^2 + rt(100, 2)) > (e. <- ellipsoidhull(d3o, ret.sq = TRUE)) 'ellipsoid' in 3 dimensions: center = ( 0.32491 101.68998 39.48045 ); squared ave.radius d^2 = 3 and shape matrix = x x 19.655 94.364 48.739 94.364 490.860 181.022 48.739 181.022 1551.980 hence, volume = 21856 > stopifnot(all.equal(e.$sqdist, + with(e., mahalanobis(d3o, center=loc, cov=cov)), + tol = 1e-13)) > d5 <- cbind(d3o, 2*abs(y)^1.5 + rt(100,3), 3*x - sqrt(abs(y))) > (e5 <- ellipsoidhull(d5, ret.sq = TRUE)) 'ellipsoid' in 5 dimensions: center = ( -0.32451 98.54780 37.33619 1973.88383 -10.81891 ); squared ave.radius d^2 = 5 and shape matrix = x x 17.8372 87.0277 8.3389 2607.9 49.117 87.0277 446.9453 -2.0502 12745.4 239.470 8.3389 -2.0502 1192.8439 2447.8 24.458 2607.9264 12745.3826 2447.8006 384472.1 7179.239 49.1172 239.4703 24.4582 7179.2 135.260 hence, volume = 191410 > tail(sort(e5$sqdist)) ## 4 values 5.00039 ... 5.0099 [1] 4.999915 5.000005 5.000010 5.000088 5.001444 5.009849 > > (e5.1e77 <- ellipsoidhull(1e77*d5)) 'ellipsoid' in 5 dimensions: center = ( -3.2451e+76 9.8548e+78 3.7336e+78 1.9739e+80 -1.0819e+78 ); squared ave.radius d^2 = 5 and shape matrix = x x 1.7837e+155 8.7028e+155 8.3389e+154 2.6079e+157 4.9117e+155 8.7028e+155 4.4695e+156 -2.0502e+154 1.2745e+158 2.3947e+156 8.3389e+154 -2.0502e+154 1.1928e+157 2.4478e+157 2.4458e+155 2.6079e+157 1.2745e+158 2.4478e+157 3.8447e+159 7.1792e+157 4.9117e+155 2.3947e+156 2.4458e+155 7.1792e+157 1.3526e+156 hence, volume = exp(898.66) > stopifnot(# proof correct scaling c^5 + all.equal(volume(e5.1e77, log=TRUE) - volume(e5, log=TRUE), + ncol(d5) * 77* log(10)) + ) > > proc.time() user system elapsed 0.117 0.034 0.210