linear_algebra_stats

Simple Linear Algebra Functions for Statistics

Description

‘broadcast’ provides some simple Linear Algebra Functions for Statistics:
cinv()
sd_lc()
ecumprob()

Usage

cinv(x)

sd_lc(X, vc, bad_rp = NaN)

ecumprob(y, sim, eps = 0)

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.

y values to estimate the cumulative probability for.

sim a matrix (or data.frame) with at least 500 columns of simulated values.
If sim is given as a dimensionless vector, it will be treated as a matrix with 1 row and length(sim) columns, and this will be noted with a message.

eps a non-negative numeric scaler smaller than 0.1, giving the cut-off value for probabilities.
Probabilities smaller than eps will be replaced with eps, and probabilities larger than 1 - eps will be replaced with 1 - eps.
Set eps = 0 to disable probability trimming.

Details

cinv()
cinv() computes the Choleski inverse of a real symmetric positive-definite square matrix.

sd_lc()
Given the linear combination X %*% b, where:

X is a matrix of multipliers/constants;
b is a vector of (correlated) random variables;
vc is the symmetric variance-covariance matrix for b;

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’).

ecumprob()
The ecumprod(y, sim) function takes a matrix (or data.frame) of simulated values sim, and for each row i (after broadcasting), estimates the cumulative distribution function of sim[i, ], and returns the cumulative probability for y[i].

In terms of statistics, it is equivalent to the following operation for each index i:
ecdf(sim[i,])(y[i])
However, ecumprob() is much faster, and supports NAs/NaNs.

In terms of linear algebra, it is equivalent to the following broadcasted operation:
rowMeans(sim <= y)
where y and sim are broadcaster arrays.
However, ecumprob() is much more memory-efficient, supports a data.frame for sim, and has statistical safety checks.

Value

For cinv():
A matrix.

For sd_lc():
A vector of standard deviations.

For ecumprob():
A vector of cumulative probabilities.
If for any observation i (after broadcasting), y[i] is NA/NaN or any of sim[i,] is NA/NaN, the result for i will be NA.
If zero-length y or sim is given, a zero-length numeric vector is returned.

References

John A. Rice (2007), Mathematical Statistics and Data Analysis (6th Edition)

Examples

library("broadcast")


# variances ====
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]  23.18594182  58.14607162  34.27339432  32.28476554  10.93475350
#>   [6]  19.24168468  36.91661748  38.92937001   0.09604826  13.55537283
#>  [11]  65.61885428   9.43946882   0.76641024  26.53762776  25.92132816
#>  [16]  26.03322670  40.15308386  17.76423215  13.85465371  19.10481584
#>  [21]   9.02798560   9.16136462  28.35501637  30.54722549  16.75492220
#>  [26]  58.14854167  36.72877757  39.83679500  22.79293876  39.77087712
#>  [31]  18.70804026  39.43864610  12.78972429  21.96297958   7.61901299
#>  [36]  43.48673728   4.51898480  38.06475711  27.08126117  47.20583632
#>  [41]  38.53069900  29.72965226   5.41031394  18.34307431  47.65382727
#>  [46]  36.18299909  42.51546063  10.13121645  19.93433892  32.72450173
#>  [51]  17.16106368   5.12765413   4.32398259   2.26509333  56.24283283
#>  [56]  18.25201245  22.64193112  27.28711613  26.98491609  35.96467012
#>  [61] 110.57966804   3.49978303  25.69103267  24.54147534   6.09460056
#>  [66]  40.36310789   4.99633511  42.04627203  36.58094067  27.91771872
#>  [71]  16.34412622  38.55372964  29.49890771  18.16415528   3.51298410
#>  [76]  57.31163256  23.22014040  26.18610584  25.92235634  11.60769305
#>  [81]  14.78182767   1.68496353  13.85877086  33.90020982   2.57156370
#>  [86]   7.74430150   1.31459335   1.71452483  22.13882479   9.09777109
#>  [91]  73.67519183  42.52566123  15.09131690  15.16787807   0.52663982
#>  [96]  46.16171218  15.59868345  38.98880181  51.06437072  24.36791662



# ecumprob() ====

sim <- rnbinom(10 * 1e4, mu = 3, size = 2) |> matrix(10, 1e4)
y <- sample(0:9)

# vector:
pnbinom(y[1], mu = 3, size = 2) # real probability
#> [1] 0.76672
ecumprob(y[1], sim[1, , drop = TRUE]) # approximation
#> [1] 0.7732

# matrix:
cbind(
  real = pnbinom(y, mu = 3, size = 2), # real probability
  approx = ecumprob(y, sim) # approximation
)
#>            real approx
#>  [1,] 0.7667200 0.7732
#>  [2,] 0.8413696 0.8449
#>  [3,] 0.1600000 0.1618
#>  [4,] 0.9294561 0.9303
#>  [5,] 0.9697669 0.9717
#>  [6,] 0.3520000 0.3520
#>  [7,] 0.6630400 0.6646
#>  [8,] 0.5248000 0.5253
#>  [9,] 0.9536426 0.9547
#> [10,] 0.8936243 0.8919

# data.frame:
cbind(
  real = pnbinom(y, mu = 3, size = 2), # real probability
  approx = ecumprob(y, as.data.frame(sim)) # approximation
)
#>            real approx
#>  [1,] 0.7667200 0.7732
#>  [2,] 0.8413696 0.8449
#>  [3,] 0.1600000 0.1618
#>  [4,] 0.9294561 0.9303
#>  [5,] 0.9697669 0.9717
#>  [6,] 0.3520000 0.3520
#>  [7,] 0.6630400 0.6646
#>  [8,] 0.5248000 0.5253
#>  [9,] 0.9536426 0.9547
#> [10,] 0.8936243 0.8919