armadillo 8.400 RC1

[conradsnicta]

No description provided.

[coatless]

@conradsnicta a thousand thanks for the work on the random number generation, e.g. scalars + wishart! (You've also just blown up the use one of my packages... But, I'll take it!)

[coatless]

@conradsnicta few quick remarks:

• Implementation wise, it matches with what I was running:

https://github.com/coatless/r-to-armadillo/blob/3ebbc49b219c3d49724c110709f4d545108f7f7f/src/distributions.cpp#L28-L55

https://github.com/conradsnicta/armadillo-code/blob/d5c0f728a0bf392ea93be10a4d34a793230b7d8c/include/armadillo_bits/op_wishrnd_meat.hpp#L134-L155

• df should always be an integer.
• Could we switch it from wishrnd(S, df) to wishrnd(df, S)? The latter is the preferred notation.
• I tested the wishart generation using my branch of RcppArmadillo and it was within the error tolerances:

# Build the C++ function
Rcpp::cppFunction("arma::mat test_rwishart(arma::mat S, int df) {
  return arma::wishrnd(S, df);
}", plugins = "cpp11", depends = "RcppArmadillo")

# Run through tests:

## Artificial
S = toeplitz((10:1)/10)

set.seed(11)
R = replicate(1000, test_rwishart(S, 20))

dim(R) 
#  10 10  1000

mR = apply(R, 1:2, mean)  
# ~= E[ Wish(S, 20) ] = 20 * S

stopifnot(all.equal(mR, 20*S, tolerance = .009))

## See Details, the variance is
Va = 20*(S^2 + tcrossprod(diag(S)))
vR = apply(R, 1:2, var)

stopifnot(all.equal(vR, Va, tolerance = 1/16))

[kthohr]

Hi @conradsnicta

I checked

• wishrnd and iwishrnd
• randg
• normcdf
• kron (for sp_mat)

and everything came out as expected.

One suggestion, however. Based on my experience, the df_simple case here

https://github.com/conradsnicta/armadillo-code/blob/unstable/include/armadillo_bits/op_wishrnd_meat.hpp#L124

is very inefficient when df is large (as tmp'tmp is costly to calculate). Example: when df = 1E05 and S = eye(3,3), the second implementation---similar to what I use in my library, but for integer df---runs ~ 100 times faster than df_simple, even with multi-threaded BLAS.

Perhaps remove the df_simple case entirely, or include a (df < some_bound) switch?

[eddelbuettel]

I ran a complete reverse depends check and just going over the results -- we may have one regression in one package that I need to look at more closely.

[eddelbuettel]

Looked at that package--it fails with the CRAN version of RcppArmadillo too so no regression here.