library("broadcast")
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] 39.1494606 48.7642509 14.0169729 5.8540772 30.5050670 16.0051528
## [7] 0.5348298 6.2267884 1.7503189 42.2695600 13.2497832 30.7276539
## [13] 22.1648623 0.9303296 36.8748101 28.0210150 7.2180810 12.6020644
## [19] 38.0135678 52.5362176 51.1606585 39.5113160 18.7700858 6.5252945
## [25] 3.2410921 6.9912552 16.2410208 20.0613612 48.1293443 24.0186162
## [31] 9.8465027 37.0471054 33.3630900 54.9061608 77.1084270 49.5606102
## [37] 7.8372873 30.6354795 57.4448232 8.4858634 20.6535663 16.4401561
## [43] 0.6167611 18.3495403 4.9164191 33.1753734 4.2067605 1.8554829
## [49] 9.5920942 3.0670750 58.4769241 26.1506793 27.4389803 44.7045079
## [55] 8.0028857 79.4254996 27.3557100 8.3145375 55.3159413 32.9990923
## [61] 74.5510744 4.1686131 27.2424800 54.3006710 4.0247830 20.7099344
## [67] 21.3683395 44.2942790 0.1471875 3.1065987 30.7900007 3.8573873
## [73] 11.7536542 6.2772705 4.5702748 16.1441767 18.8832618 51.1590735
## [79] 18.3127694 15.8591069 13.3458475 60.1090373 4.4611244 7.3424929
## [85] 10.7595152 61.8141969 51.8672997 20.8724230 10.6810167 8.1468471
## [91] 28.7656320 39.6949603 9.1058719 10.7463119 0.0586897 32.7616679
## [97] 34.1716676 80.2903029 5.6503084 9.0225351linear_algebra_stats
Simple Linear Algebra Functions for Statistics
Description
‘broadcast’ provides some simple Linear Algebra Functions for Statistics:
cinv();
sd_lc().
Usage
cinv(x)
sd_lc(X, vc, bad_rp = NaN)
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. |
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’).
Value
For cinv():
A matrix.
For sd_lc():
A vector of standard deviations.
References
John A. Rice (2007), Mathematical Statistics and Data Analysis (6th Edition)
See Also
chol, chol2inv