library("broadcast")
# variances ====
vc <- datasets::ability.cov$cov
X <- matrix(rnorm(100), 100, ncol(vc))
solve(vc)
#> general picture blocks maze reading
#> general 0.082259001 -0.0312406436 -7.750932e-03 -0.013309494 -2.061705e-02
#> picture -0.031240644 0.2369906996 -2.484938e-02 0.017844845 8.603286e-04
#> blocks -0.007750932 -0.0248493822 1.344272e-02 -0.012544830 -3.802671e-05
#> maze -0.013309494 0.0178448450 -1.254483e-02 0.101625400 5.508423e-03
#> reading -0.020617049 0.0008603286 -3.802671e-05 0.005508423 5.713620e-02
#> vocab -0.002420800 0.0019394999 -1.157864e-03 -0.002857265 -2.406969e-02
#> vocab
#> general -0.002420800
#> picture 0.001939500
#> blocks -0.001157864
#> maze -0.002857265
#> reading -0.024069692
#> vocab 0.020323179
cinv(vc) # faster than `solve()`, but only works on positive definite matrices
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.082259001 -0.0312406436 -7.750932e-03 -0.013309494 -2.061705e-02
#> [2,] -0.031240644 0.2369906996 -2.484938e-02 0.017844845 8.603286e-04
#> [3,] -0.007750932 -0.0248493822 1.344272e-02 -0.012544830 -3.802671e-05
#> [4,] -0.013309494 0.0178448450 -1.254483e-02 0.101625400 5.508423e-03
#> [5,] -0.020617049 0.0008603286 -3.802671e-05 0.005508423 5.713620e-02
#> [6,] -0.002420800 0.0019394999 -1.157864e-03 -0.002857265 -2.406969e-02
#> [,6]
#> [1,] -0.002420800
#> [2,] 0.001939500
#> [3,] -0.001157864
#> [4,] -0.002857265
#> [5,] -0.024069692
#> [6,] 0.020323179
all(round(solve(vc), 6) == round(cinv(vc), 6)) # they're the same
#> [1] TRUE
sd_lc(X, vc)
#> [1] 8.17822152 8.52950790 14.61061866 13.74182009 20.08895593 52.74320344
#> [7] 13.18236441 32.13640467 62.07103664 26.10621542 27.75234330 8.87233237
#> [13] 5.99264222 27.46148581 23.36506122 72.63498851 28.36194356 40.85363451
#> [19] 2.89106771 70.36078557 18.44520520 40.27459994 2.58513738 12.15714633
#> [25] 42.40624802 3.14036425 6.78902216 6.63237839 48.71984088 71.04784554
#> [31] 69.62394097 25.05717117 19.83321567 51.30543443 24.09174345 15.03181068
#> [37] 8.59054261 15.59001194 24.11932712 19.06114013 8.33432389 24.87198679
#> [43] 0.07269481 56.85713632 19.31126104 31.36083876 1.69428990 18.89291174
#> [49] 1.19857516 7.52815451 58.52979580 22.54981646 18.94014026 35.28797336
#> [55] 32.92857221 12.85891059 31.41966510 76.73736264 39.75320360 22.24164083
#> [61] 86.66226633 40.52532750 9.62719848 52.57748298 14.25724314 9.62252908
#> [67] 9.57276083 13.35340190 10.92950851 5.57750898 50.01780201 9.92045810
#> [73] 0.51821409 15.19833344 77.31974983 11.65728049 46.32635869 22.16827101
#> [79] 2.42364970 5.79950255 13.13485718 16.67429518 0.03232121 1.21111082
#> [85] 48.43809634 10.82356982 0.39740737 38.18113426 76.33010790 2.72615563
#> [91] 24.73730057 19.73174845 46.71974062 12.78622339 0.75745292 0.01928867
#> [97] 22.06790892 9.56908640 40.29155399 14.20337682
# ecumprob() ====
sim <- rnbinom(10 * 5000L, mu = 3, size = 2) |> matrix(10, 5000)
y <- sample(0:9)
# vector:
pnbinom(y[1], mu = 3, size = 2) # real probability
#> [1] 0.9697669
ecumprob(y[1], sim[1, , drop = TRUE]) # approximation
#> [1] 0.9694
# matrix:
pnbinom(y, mu = 3, size = 2) # real probability
#> [1] 0.9697669 0.8936243 0.5248000 0.1600000 0.6630400 0.3520000 0.9536426
#> [8] 0.7667200 0.9294561 0.8413696
ecumprob(y, sim) # approximation
#> [1] 0.9694 0.8922 0.5254 0.1526 0.6662 0.3622 0.9526 0.7704 0.9318 0.8372
# data.frame:
pnbinom(y, mu = 3, size = 2) # real probability
#> [1] 0.9697669 0.8936243 0.5248000 0.1600000 0.6630400 0.3520000 0.9536426
#> [8] 0.7667200 0.9294561 0.8413696
ecumprob(y, as.data.frame(sim)) # approximation
#> [1] 0.9694 0.8922 0.5254 0.1526 0.6662 0.3622 0.9526 0.7704 0.9318 0.8372linear_algebra_stats
Simple Linear Algebra Functions for Statistics
Description
‘broadcast’ provides some simple Linear Algebra Functions for Statistics:
cinv()
sd_lc()
ecumprob()
Usage
cinv(x)
sd_lc(X, vc, bad_rp = NaN)
ecumprob(y, sim, eps = 0)
Arguments
x
|
a real symmetric positive-definite square matrix. |
X
|
a numeric (or logical) matrix of multipliers/constants |
vc
|
the variance-covariance matrix for the (correlated) random variables. |
bad_rp
|
if vc is not a Positive (semi-) Definite matrix, give here the value to replace bad standard deviations with. |
y
|
values to estimate the cumulative probability for. |
sim
|
a matrix (or data.frame) with at least 500 columns of simulated values. If sim is given as a dimensionless vector, it will be treated as a matrix with 1 row and length(sim) columns, and this will be noted with a message.
|
eps
|
a non-negative numeric scaler smaller than 0.1, giving the cut-off value for probabilities. Probabilities smaller than eps will be replaced with eps, and probabilities larger than 1 - eps will be replaced with 1 - eps. Set eps = 0 to disable probability trimming.
|
Details
cinv()
cinv() computes the Choleski inverse of a real symmetric positive-definite square matrix.
sd_lc()
Given the linear combination X %*% b, where:
-
Xis a matrix of multipliers/constants; -
bis a vector of (correlated) random variables; -
vcis the symmetric variance-covariance matrix forb;
sd_lc(X, vc) computes the standard deviations for the linear combination X %*% b, without making needless copies.
sd_lc(X, vc) will use much less memory than a base ‘R’ approach.
sd_lc(X, vc) will usually be faster than a base ‘R’ approach (depending on the Linear Algebra Library used for base ‘R’).
ecumprob()
The ecumprod(y, sim) function takes a matrix (or data.frame) of simulated values sim, and for each row i (after broadcasting), estimates the cumulative distribution function of sim[i, ], and returns the cumulative probability for y[i].
In terms of statistics, it is equivalent to the following operation for each index i:
ecdf(sim[i,])(y[i])
However, ecumprob() is much faster, and supports NAs/NaNs.
In terms of linear algebra, it is equivalent to the following broadcasted operation:
rowMeans(sim <= y)
where y and sim are broadcaster arrays.
However, ecumprob() is much more memory-efficient, supports a data.frame for sim, and has statistical safety checks.
Value
For cinv():
A matrix.
For sd_lc():
A vector of standard deviations.
For ecumprob():
A vector of cumulative probabilities.
If for any observation i (after broadcasting,) y[i] is NA/NaN or any of sim[i,] is NA/NaN, the result for i will be NA.
If zero-length y or sim is given, a zero-length numeric vector is returned.
References
John A. Rice (2007), Mathematical Statistics and Data Analysis (6th Edition)
See Also
chol, chol2inv