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.323605 24.021107 8.536008 4.276802 8.959474 48.316251 5.380940
## [8] 12.695256 79.514488 7.805131 22.651060 4.609553 32.445886 16.786096
## [15] 21.562626 23.694818 46.299976 12.916080 14.052731 30.469274 5.964171
## [22] 14.863878 33.196687 2.456307 56.656684 54.263042 16.744010 9.224998
## [29] 18.988468 24.434713 22.452238 28.106863 40.801421 23.688360 7.502245
## [36] 18.852240 16.250687 18.861764 5.775901 16.941047 32.300444 23.402769
## [43] 1.288353 40.264223 42.580749 16.810504 22.682131 1.158418 63.671497
## [50] 42.142986 15.534527 4.748445 48.704755 62.661159 24.930954 44.405448
## [57] 34.065715 60.227810 23.316176 11.874941 39.020084 9.077244 6.111850
## [64] 9.882684 6.976850 28.072847 18.987378 20.547105 8.812027 30.753153
## [71] 26.815797 49.740254 32.723872 15.784515 34.847668 35.781363 20.132517
## [78] 12.007347 25.587716 44.307504 4.361693 5.060369 31.125893 39.533376
## [85] 52.845217 16.382416 4.997381 19.997515 12.281055 55.792841 4.571092
## [92] 11.003290 22.530021 20.119462 17.051960 36.282227 14.278220 40.793655
## [99] 59.428830 34.859806linear_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) may possibly, but not necessarily, 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