I’ve written a couple packages for simulating from some distributions (Wishart-related distributions in CholWishart and matrix-variate distributions in matrixdist) and sometimes when a new function has been written or has been refactored, you need some ways to verify it is giving the answers you expect. There is an entire literature on this in general, but I am going to discuss one handy trick I sometimes use if there isn’t much theory to rely on or if you just want some quick heuristic to show you are in the right ballpark.
Sometimes there are known theoretical distributions of test statistics based on your distribution, in which case those can be used (if they don’t rely too heavily on asymptotics and can actually distinguish from something wrong).
The way I implemented rCholWishart and rInvWishart, they should give, with the same seeds, equivalent answers to rWishart if an identity covariance matrix is used, so this becomes fairly easy: give them the same seeds and see if they are what is expected.
library('CholWishart')
set.seed(20190109)
A <- rWishart(1, 5, diag(3))[,,1]
set.seed(20190109)
B <- rInvWishart(1, 5, diag(3))[,,1]
A %*% B## [,1] [,2] [,3]
## [1,] 1.000000e+00 0.000000e+00 4.440892e-16
## [2,] 0.000000e+00 1.000000e+00 1.110223e-16
## [3,] 4.440892e-16 -5.551115e-17 1.000000e+00Incidentally, the package bench is useful if you want to benchmark two versions of the same function which should output the exact same thing.
However, some of my functions are re-implementations of functions in other packages that do not use the results of the random number generator in the same way, so, while I trust (or should I?) their results, mine are not exactly the same. Or perhaps I am trying an efficient re-implementation of a slow function and need to verify the results are as they seem.
Consider this in the case of the pseudo Wishart. Let \(X_i \sim MVN_p(\mathbf{0}, \Sigma)\) for \(i = 1, 2, \ldots, n\) be \(n\) i.i.d. draws from a \(p\)-dimensional multivariate normal distribution with mean \(\mathbf{0}\) and covariance matrix \(\Sigma\). Then if \(n > p -1\), the \(\sum X_i X_i^T\) is distributed as a Wishart random matrix with \(n\) degrees of freedom. If \(n \leq p-1\), this is a singular distribution known as the pseudo Wishart distribution.
There is a function in another package, rWishart, which handles it. I wrote a slow version myself for my own package (v0.9.4) and a fast version for the latest (still in development) version (v1.0.0). We may want to verify that things are what they seem. It is straightforward to simulate from this distribution: simulate from a multivariate normal (try MASS::mvrnorm), multiply them together appropriately, and there is your random matrix. This can serve as a “ground truth” if we trust MASS and our multiplicative abilities. Then we can compare to another package (presumed to be correct) and then to the function that is going to be in this package. We should not be surprised if the naive “ground truth” is faster than the version in some package, since this version does not include input validation or error checking, which a function for general use should.
What I like to do as a first sanity check is to simulate a modest number (say, 1000) and make sure the figures are roughly comparable. Note, by the way, I start by simulating a covariance matrix to use as input to the functions - I don’t want to use an identity matrix because we want to know that the things are factoring the covariance matrices appropriately. So here we are going to simulate 1000 draws from a pseudo-Wishart with 3 degrees of freedom and dimension 5 with a covariance matrix \(A\). If we wanted to, we could profile these results and compare timings as well.
A <- rWishart(1, 10, diag(5))[,,1]
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 9.2262807 -6.200813 -0.8685732 3.006720 -7.395744
## [2,] -6.2008125 8.927408 0.6052070 2.317195 7.305159
## [3,] -0.8685732 0.605207 8.0806719 -3.054297 -0.457451
## [4,] 3.0067197 2.317195 -3.0542970 11.574879 -1.452411
## [5,] -7.3957440 7.305159 -0.4574510 -1.452411 9.749766# There are better ways of doing B, but I want a quick-and-dirty
# that is verifiably correct
B <- MASS::mvrnorm(3000, rep(0,5), A)
Blist <- list()
for (i in 1:1000) Blist[[i]] <- crossprod(B[(3*i-2):(3*i),])
C <- rWishart::rPsuedoWishart(1000, 3, A)
D <- CholWishart::rPseudoWishart(1000, 3, A)The first thing I would do is compare the sum of the elements - these should be roughly the same. Here there is a theoretical result - the sum should be distributed as a Wishart with 3000 degrees of freedom.
Reduce(`+`, Blist)## [,1] [,2] [,3] [,4] [,5]
## [1,] 27565.321 -18651.829 -2138.801 8031.953 -21527.564
## [2,] -18651.829 26757.781 1388.668 7527.585 21613.654
## [3,] -2138.801 1388.668 24792.339 -9753.315 -2073.009
## [4,] 8031.953 7527.585 -9753.315 34720.093 -3137.397
## [5,] -21527.564 21613.654 -2073.009 -3137.397 28267.008rowSums(x = C, na.rm = FALSE, dims = 2L)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1490828.55 -469207.7 -52086.42 265758.24 -282072.87
## [2,] -469207.75 455437.8 -145829.39 135607.77 107591.31
## [3,] -52086.42 -145829.4 285730.29 -128320.17 -15181.54
## [4,] 265758.24 135607.8 -128320.17 216503.43 -42798.66
## [5,] -282072.87 107591.3 -15181.54 -42798.66 62265.16rowSums(D, FALSE, 2L)## [,1] [,2] [,3] [,4] [,5]
## [1,] 27556.89 -18049.846 -2081.950 9223.790 -22020.422
## [2,] -18049.85 26323.893 1137.398 7193.153 21898.722
## [3,] -2081.95 1137.398 23471.798 -8921.188 -1738.466
## [4,] 9223.79 7193.153 -8921.188 35459.581 -4111.849
## [5,] -22020.42 21898.722 -1738.466 -4111.849 29202.886The expected value of these is \(3000 A\) and you can find the variance of individual terms if you really want, but this is perhaps better as a rough sanity check during development. Here it can be seen that the three methods have roughly the same average result. Not shown here, but you may also which to compare, say, box plots of the distributions of the 25 elements and repeat this with some other covariance matrices as input (say, an \(AR(1)\) matrix). This exercise may seem a little trivial, but, for instance, the inverse Wishart has a few different definitions of the degrees of freedom parameter. If you wish to substitute one function for another (that is faster), it is important to verify that they are parameterized in the same way or that the reparameterization is the same.
Laplacelist <- list()
for (i in 1:1000) Laplacelist[[i]] <- LaplacesDemon::rinvwishart(10, A)
MyInvWishart <- CholWishart::rInvWishart(1000, 10, A)
Reduce(`+`, Laplacelist)## [,1] [,2] [,3] [,4] [,5]
## [1,] 2305.8541 -1571.6677 -217.7485 648.9478 -1852.5833
## [2,] -1571.6677 2247.8176 203.7237 646.2970 1824.6910
## [3,] -217.7485 203.7237 2202.1261 -796.8863 -163.5460
## [4,] 648.9478 646.2970 -796.8863 2867.5037 -254.8008
## [5,] -1852.5833 1824.6910 -163.5460 -254.8008 2449.3090rowSums(MyInvWishart, FALSE, 2L)## [,1] [,2] [,3] [,4] [,5]
## [1,] 2253.3570 -1537.0797 -291.71139 767.0758 -1820.76919
## [2,] -1537.0797 2250.1777 181.11747 614.5197 1840.93160
## [3,] -291.7114 181.1175 1969.53377 -776.3872 -41.93022
## [4,] 767.0758 614.5197 -776.38720 3017.8656 -347.40198
## [5,] -1820.7692 1840.9316 -41.93022 -347.4020 2429.55687This suggests, then, that the results functions have the same parameterization for the degrees of freedom.