Issues sampling from Wishart Distribution

Hi: I was wondering if there were changes made to the functions, as it does not seem to work for me as expected. More specifically, if taking the average of the diagonals everything appears OK, but each diagonal should be 20. This is note the case for the arma implementation, and the values can be very far from the expected (which also seems influenced by the correlation sizes for some reason).

R solution

samples = 10000
S = toeplitz((10:1)/10)

R_wish <- rWishart(samples, 20, S)
R_list <- list()
for(i in 1:samples){
  R_list[[i]] <- diag(R_wish[,,i])
}

colMeans(do.call(rbind, R_list) )
20.05181 19.96682 19.93232 19.91842 19.89735 19.86378 19.87150 19.90565 19.91630 19.92229

# now take the average of the diagonals as
mean(colMeans(do.call(rbind, R_list) ))

19.92462

Rcpp solution

Rcpp::cppFunction("arma::mat test_rwishart(arma::mat S, int df) {
  return arma::wishrnd(S, df);
                  }", plugins = "cpp11", depends = "RcppArmadillo")


Rcpp_wish = replicate(samples, test_rwishart(S, 20))


Rcpp_list <- list()
for(i in 1:samples){
 Rcpp_list[[i]] <- diag(Rcpp_wish[,,i])
  
}
colMeans( do.call(rbind, Rcpp_list))

# these should all be 20
77.111685 21.990086 20.196240 18.345458 16.408652 14.297822 11.994180  9.526644  6.773458  3.668038

# now take the average of the diagonals as

mean(colMeans( do.call(rbind, Rcpp_list)))

20.03123

Thanks !