R Under development (unstable) (2024-01-12 r85803 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. > options(warn=1) > require(compositions) Loading required package: compositions Welcome to compositions, a package for compositional data analysis. Find an intro with "? compositions" Attaching package: 'compositions' The following objects are masked from 'package:stats': anova, cor, cov, dist, var The following object is masked from 'package:graphics': segments The following objects are masked from 'package:base': %*%, norm, scale, scale.default > > res <- rAitchison(100,theta=c(1,2,3),sigma=ilrvar2clr(diag(c(0.1,2)))) > res2<- rnorm.acomp(100,acomp(c(1,2,3)),ilrvar2clr(diag(c(0.1,2)))) > plot(res) > plot(res2,add=TRUE,col="red") > > dr = dAitchison(res,theta=c(1,2,3),sigma=ilrvar2clr(diag(c(1,2)))) > print(dr) [1] 0.103599321 0.097915213 0.100882901 0.092219836 0.013763681 0.065324109 [7] 0.017851278 0.093417797 0.044606617 0.085452660 0.080984272 0.021076987 [13] 0.047004299 0.004107832 0.045285416 0.041288857 0.039527487 0.023656865 [19] 0.010923871 0.026988504 0.013972007 0.053561572 0.046701581 0.099720605 [25] 0.091230601 0.088775806 0.089372430 0.079629281 0.073764734 0.098887068 [31] 0.076915249 0.094159381 0.059215042 0.058570513 0.090326158 0.059157553 [37] 0.053187646 0.032238059 0.088034866 0.019570246 0.093042856 0.101380414 [43] 0.037820140 0.030865802 0.097298433 0.060484320 0.099065495 0.048123182 [49] 0.029520329 0.044908779 0.089059500 0.065538870 0.033319486 0.100613026 [55] 0.057335898 0.028197555 0.098230458 0.095073927 0.091919004 0.084154092 [61] 0.079271289 0.073940073 0.059999352 0.028661052 0.052196255 0.072074034 [67] 0.076357507 0.001920746 0.087571521 0.094355811 0.089341769 0.043275221 [73] 0.104890581 0.007716587 0.070136307 0.026420461 0.002194903 0.058214978 [79] 0.074121535 0.065009337 0.095325638 0.023952403 0.074445187 0.095676738 [85] 0.038646484 0.064735819 0.098052295 0.022728492 0.100291561 0.074310258 [91] 0.086373240 0.067379271 0.016080817 0.005196246 0.072355828 0.081826667 [97] 0.088744982 0.032267845 0.037760807 0.060735245 > > > > erg<-AitchisonDistributionIntegrals(c(-1,3,-2),ilrvar2clr(-diag(c(1,2))),grid=60) > print(erg) $theta [1] -1 3 -2 $beta [,1] [,2] [,3] [1,] -0.8333333 0.1666667 0.6666667 [2,] 0.1666667 -0.8333333 0.6666667 [3,] 0.6666667 0.6666667 -1.3333333 $alpha [1] 0 $mu [1] "0.1033797" "0.7638782" "0.1327421" attr(,"class") [1] "acomp" $sigma [,1] [,2] [,3] [1,] 0.29166667 -0.20833333 -0.08333333 [2,] -0.20833333 0.29166667 -0.08333333 [3,] -0.08333333 -0.08333333 0.16666667 $expKappa [1] 126.4483 $loggxMean [1] -1.595297 $clrMean [1] -0.7491284 1.2492295 -0.5001012 $clrVar [,1] [,2] [,3] [1,] 0.28804158 -0.20476955 -0.08327203 [2,] -0.20476955 0.28841772 -0.08364817 [3,] -0.08327203 -0.08364817 0.16692020 > > (myvar<-with(erg, -1/2*ilrvar2clr(solve(clrvar2ilr(beta))))) 1 2 3 1 0.29166667 -0.20833333 -0.08333333 2 -0.20833333 0.29166667 -0.08333333 3 -0.08333333 -0.08333333 0.16666667 > (mymean<-with(erg,myvar%*%theta)) [,1] 1 -0.75 2 1.25 3 -0.50 > > with(erg,myvar-clrVar) 1 2 3 1 3.625083e-03 -0.0035637819 -6.130155e-05 2 -3.563782e-03 0.0032489448 3.148371e-04 3 -6.130155e-05 0.0003148371 -2.535355e-04 > with(erg,mymean-clrMean) [,1] 1 -0.0008716054 2 0.0007704529 3 0.0001011525 > > > AitchisonDistributionIntegrals function (theta = alpha + sigma %*% clr(mu), beta = -1/2 * gsi.svdinverse(sigma), alpha = mean(theta), mu = clrInv(c(sigma %*% (theta - alpha))), sigma = -1/2 * gsi.svdinverse(beta), grid = 30, mode = 3) { if (!xor(missing(theta), (missing(mu) || missing(alpha)))) stop("Specify either theta or mu and alpha") if (!xor(missing(beta), missing(sigma))) stop("Specify either beta or sigma") D <- length(theta) if (nrow(beta) == D - 1) beta <- ilrvar2clr(beta) if (any(abs(beta - t(beta)) > 1e-10)) { warning("AitchisonDistributionIntegrals: beta was not symmetric 1") print(beta) } gsiInt(dim(beta), 2) stopifnot(length(dim(beta)) == 2, ncol(beta) == D, nrow(beta) == D) erg <- .C(gsiAitchisonDistributionIntegral, D = gsiInt(D, 1), grid = gsiInt(grid, 1), mode = gsiInt(mode, 1), theta = gsiDouble(theta, D), beta = gsiDouble(beta, D * D), expKappa = numeric(1), loggxMean = numeric(1), clrMean = numeric(D), clrVar = numeric(D * D)) erg$beta <- matrix(erg$beta, nrow = D) erg$SqIntegral <- matrix(erg$clrVar, nrow = D, ncol = D) erg$alpha = alpha erg$mu = mu erg$sigma = sigma erg$clrSqExpectation <- matrix(erg$clrVar, nrow = D, ncol = D) dim(erg$clrVar) <- c(D, D) return(erg[c("theta", "beta", "alpha", "mu", "sigma", if (mode >= 0) c("expKappa", "loggxMean") else c(), if (mode >= 1) "clrMean" else c(), if (mode == 2) "clrSqExpectation" else if (mode >= 3) "clrVar" else c())]) } > > > res <- rAitchison(100,theta=c(0.5,1,3),sigma=ilrvar2clr(diag(c(0.1,2)))) > plot(res) > AitStats = AitchisonDistributionIntegrals(theta=c(0.5,1,3),sigma=ilrvar2clr(diag(c(0.1,2)))) > plot(clrInv(AitStats$clrMean), add=TRUE, col=2, pch=19) > print(AitStats) $theta [1] 0.5 1.0 3.0 $beta [,1] [,2] [,3] [1,] -2.54166667 2.45833333 0.08333333 [2,] 2.45833333 -2.54166667 0.08333333 [3,] 0.08333333 0.08333333 -0.16666667 $alpha [1] 1.5 $mu [1] "0.01059915" "0.01114258" "0.97825827" attr(,"class") [1] "acomp" $sigma 1 2 3 1 0.3833333 0.2833333 -0.6666667 2 0.2833333 0.3833333 -0.6666667 3 -0.6666667 -0.6666667 1.3333333 $expKappa [1] 0.08158924 $loggxMean [1] -1.346044 $clrMean [1] -0.3959256 -0.3503244 0.7462500 $clrVar [,1] [,2] [,3] [1,] 0.1345549 0.0430842 -0.1776391 [2,] 0.0430842 0.1340438 -0.1771280 [3,] -0.1776391 -0.1771280 0.3547671 > > proc.time() user system elapsed 0.45 0.04 0.48