R Under development (unstable) (2024-01-07 r85787 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. > > 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.080902627 0.085382639 0.094304883 0.012804406 0.053435804 0.074763754 [7] 0.036116138 0.090484115 0.047089602 0.078531764 0.081834546 0.039073909 [13] 0.069214969 0.086782019 0.085999448 0.024209419 0.082757704 0.048025437 [19] 0.082812683 0.102463961 0.037684483 0.063250767 0.092208066 0.094525700 [25] 0.077762013 0.093319031 0.092518858 0.030658402 0.045249704 0.052992998 [31] 0.096462593 0.066294447 0.048087293 0.039723290 0.061359951 0.080547535 [37] 0.093146276 0.045905326 0.097517852 0.097392217 0.083939758 0.098066322 [43] 0.092225638 0.018415917 0.023553366 0.064349067 0.023053541 0.040850091 [49] 0.052206175 0.042186902 0.041780630 0.053081460 0.095485779 0.001739205 [55] 0.052533965 0.040442444 0.003291406 0.053855692 0.013214769 0.076476452 [61] 0.061673325 0.102774141 0.069495419 0.087471191 0.065616175 0.052333115 [67] 0.028947179 0.090520138 0.047636873 0.047827635 0.079639523 0.044663898 [73] 0.073580214 0.085371023 0.104506861 0.018061836 0.077640281 0.075789830 [79] 0.098495537 0.065331107 0.101419548 0.078909023 0.065595606 0.085471006 [85] 0.067522217 0.099445728 0.092335247 0.033298394 0.097988313 0.040062411 [91] 0.077379667 0.042699103 0.104405384 0.088878070 0.068678563 0.069285975 [97] 0.043357534 0.075556848 0.059008019 0.040673277 > > > > 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.42 0.09 0.50