R Under development (unstable) (2025-08-19 r88650 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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.062918249 0.102087762 0.089940491 0.010863076 0.103870904 0.010110156 [7] 0.075141939 0.051053340 0.091090068 0.096744171 0.026571288 0.084647965 [13] 0.063346586 0.099723258 0.104636825 0.053541534 0.045990325 0.064027410 [19] 0.101772157 0.095226531 0.100196311 0.096251133 0.093932532 0.059949911 [25] 0.089121963 0.095568898 0.068875757 0.072304315 0.020411088 0.044054908 [31] 0.064104197 0.078171153 0.097523277 0.078253649 0.041480477 0.077932778 [37] 0.102403744 0.010301204 0.072942861 0.012286379 0.014068305 0.012636825 [43] 0.008474344 0.099516609 0.070402743 0.076872375 0.029541095 0.077903421 [49] 0.053644842 0.076404154 0.029547065 0.103585484 0.104288275 0.101076572 [55] 0.013282800 0.086607618 0.086890960 0.052856629 0.054488164 0.096360620 [61] 0.100388865 0.094502829 0.055583020 0.094582207 0.021856848 0.104881532 [67] 0.085115778 0.062749953 0.021486124 0.004878872 0.025828127 0.080286736 [73] 0.082874272 0.071979671 0.046856262 0.101299109 0.082189674 0.066831354 [79] 0.041778541 0.011204827 0.035039741 0.087045454 0.036759545 0.090984975 [85] 0.088223164 0.035256364 0.098769666 0.018312986 0.098696876 0.033026529 [91] 0.023187359 0.095480063 0.040360013 0.034527461 0.050769964 0.059553555 [97] 0.100145578 0.057884142 0.065867961 0.089149628 > > > > 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.68 0.15 0.82