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] 19.3753824 6.1507576 15.8746647 9.0430574 5.1390576 35.7575040
## [7] 17.6120319 0.5259112 2.2080528 40.2762955 11.7751464 34.7859871
## [13] 10.0858898 30.8060634 4.8281304 24.9319258 11.3993408 2.5318956
## [19] 46.6877859 28.3051241 13.8141694 27.5798488 66.7363468 45.4445903
## [25] 1.8670598 8.2352032 6.9459885 68.5198946 4.6736461 40.5609897
## [31] 52.6092793 0.8169092 10.2980189 29.7404485 56.7786546 58.5301492
## [37] 64.0120315 36.3352199 78.7263934 30.9759612 46.0931623 13.9811387
## [43] 52.1223739 17.5603817 43.1527029 24.7441478 18.3533967 16.2241645
## [49] 27.2878819 22.1219271 1.6625682 18.5339956 16.8404374 7.9373964
## [55] 43.0868674 12.6066415 37.8479873 30.6471430 3.4791505 37.6073421
## [61] 36.3911462 22.9167872 9.2736564 4.3509642 18.9158370 36.3820304
## [67] 1.7776836 13.3255258 6.8026834 10.9968654 3.5852535 85.2307616
## [73] 4.6954430 2.1895021 39.8760421 4.8920657 2.8613296 39.7972239
## [79] 36.4444679 67.1729134 15.5724761 30.4682563 29.9877194 65.0069704
## [85] 25.6310360 42.2253751 2.9518456 0.1591432 35.3947077 3.3169854
## [91] 9.3326001 17.5836413 17.3617417 23.7843456 10.2158799 16.2675801
## [97] 14.2412838 17.4214347 13.1882574 45.6662292linear_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