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] 3.1173396 18.4834052 11.0998195 34.2289288 14.1711067 20.5281797
## [7] 20.8883847 72.9521942 34.6393005 0.7284870 22.0046437 3.9454853
## [13] 3.0363658 27.1181195 16.0294652 22.6444920 42.5850437 48.9989440
## [19] 46.6109678 16.5316134 21.0861319 0.9405498 3.8020157 13.9734753
## [25] 15.4668460 33.3142207 27.2569932 28.8245209 39.5672922 29.4038231
## [31] 22.4848824 0.6125846 25.7240940 0.8300127 47.5212929 8.3672136
## [37] 17.1537155 28.5032539 26.1939141 4.7259071 28.7680155 16.1492060
## [43] 29.2691760 19.0831302 20.5312729 16.6629599 74.2507242 2.5332206
## [49] 37.3654003 72.7701207 49.9818692 20.5269574 34.9370033 43.6870127
## [55] 17.7562481 32.8163376 30.3480722 13.9911386 18.2328310 29.3917502
## [61] 31.3707319 3.3925721 5.8079890 22.1924365 87.9081033 3.8885877
## [67] 26.9565554 13.8363516 10.7916792 78.8055970 23.0542098 22.2042839
## [73] 41.9950731 5.6494288 4.7270681 7.0021605 33.4468620 8.0439060
## [79] 9.5409394 3.7688937 42.8111119 7.9118229 20.6491195 51.5434755
## [85] 27.1147114 10.8778116 5.9575351 32.2885749 32.1117972 3.1597407
## [91] 8.7573384 0.6681906 20.2307886 40.4665452 26.6282723 15.8375060
## [97] 44.2579520 9.5331960 28.5530099 32.0349369linear_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