R Under development (unstable) (2024-01-24 r85824 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.0810577335 0.0463362373 0.0617184481 0.1029478823 0.0498786828 [6] 0.0763690061 0.0714293095 0.0407974834 0.0750887033 0.0414056796 [11] 0.0613057136 0.0609635244 0.0751299138 0.0775576239 0.0201869613 [16] 0.1045156471 0.0778979166 0.0125242931 0.0551124615 0.0142219803 [21] 0.0473308595 0.0722781729 0.0364207865 0.0200310677 0.0808121354 [26] 0.0820583990 0.0795724315 0.0943346744 0.0590056309 0.1014974641 [31] 0.0627966238 0.0663133022 0.0716962321 0.0679113458 0.0287105013 [36] 0.1010447729 0.0948963330 0.0965295001 0.0130715853 0.0007039244 [41] 0.0829899094 0.1018361150 0.0850122101 0.0394696481 0.1017632942 [46] 0.0363069559 0.0754704420 0.0204075900 0.0843141737 0.0850728755 [51] 0.0504894225 0.0881230702 0.0885630540 0.0074379975 0.0935661486 [56] 0.0539078705 0.1008033551 0.0348138550 0.0322165660 0.0261531375 [61] 0.0434768485 0.0205838265 0.0241368397 0.0536268678 0.0940173227 [66] 0.0841152964 0.0404965232 0.0800213830 0.0673065087 0.0469811965 [71] 0.0656034833 0.0869137600 0.0728407784 0.0882497200 0.0481913414 [76] 0.0889148933 0.0780688615 0.0313257175 0.0606636664 0.0163188712 [81] 0.0674000485 0.0868519725 0.1002166272 0.0664890596 0.0408098232 [86] 0.1011333586 0.0652495984 0.0251359948 0.0504704067 0.0802630837 [91] 0.0855346505 0.0931968498 0.0552313276 0.0936805463 0.0088934570 [96] 0.0039747995 0.0945416697 0.0185972139 0.1015362751 0.0373950982 > > > > 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.48 0.01 0.50