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] 33.8790637 6.8419010 26.7972808 5.8714155 15.2963978 10.8256759
## [7] 17.3693665 4.2848130 20.2561527 8.4720627 41.7539987 12.7084000
## [13] 9.1997388 38.1819430 15.4644306 58.4567912 12.0821657 10.1718918
## [19] 34.5241531 39.2248695 36.2840080 42.1406248 31.6204209 4.3328182
## [25] 45.9386485 1.7025666 31.2368314 9.3482983 46.3769021 21.6163561
## [31] 9.3825170 34.4596798 8.0546023 22.1945612 14.5647937 25.5694760
## [37] 12.3332683 25.8621999 34.0441167 13.6956938 15.6243097 12.0947728
## [43] 49.6546713 4.6997393 21.9258407 2.0811662 7.4977602 43.0399231
## [49] 41.8753701 16.4755624 31.5646424 20.3113878 9.3862704 15.8894822
## [55] 23.7455341 30.4879190 10.9467654 37.9143961 17.3002746 30.3760069
## [61] 48.4928635 57.9889245 23.4071506 1.0990211 0.3867396 16.7909548
## [67] 10.9796887 37.1234765 5.8182849 0.3807780 48.4331823 31.5231952
## [73] 2.1238523 4.8507317 9.2471567 0.4488260 30.0057378 33.0544371
## [79] 40.0907926 33.7224190 5.4169435 4.6387404 36.6025119 70.8052904
## [85] 37.6979894 38.2259857 6.8312982 80.9123821 7.8221809 47.4507653
## [91] 1.6597103 15.0453860 23.5198331 47.2358615 13.0596747 7.4966369
## [97] 17.4149147 6.9269631 6.1549509 1.0878087linear_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