Wishart Distribution

> W <- matrix(c(2, -0.3, -0.3, 4), 2, 2)
> S <- matrix(c(1, -0.3, -0.3, 1), 2, 2)
> WishartPDF <- function(W, nu, S) {
+     N <- dim(W)[1]
+     ARGS <- (nu - N + 1):nu
+     K <- 2^(nu * N/2) * pi^(N * (N - 1)/4) * prod(gamma(ARGS/2))
+     f <- 1/K * (det(S))^(-nu/2) * (det(W))^((nu - N - 1)/2) *
+         exp(-1/2 * sum(diag(inv(S) %*% W)))
+     return(f)
+ }
> WishartPDF(W, nu, S)
[1] 0.0001798426

> library(MCMCpack)
> dwish(W, nu, S)
[1] 0.0001798426
> rwish(nu, S)
          [,1]      [,2]
[1,] 13.223421 -4.791452
[2,] -4.791452  5.423237

> library(mnormt)
> s <- c(1, 1)
> r <- 0.3
> nu <- 10
> Simulations <- 10000
> Sigma <- diag(s) %*% matrix(c(1, r, r, 1), nrow = 2) %*% diag(s)
> W_xx <- vector()
> W_xy <- vector()
> W_yy <- vector()
> Vec_W <- matrix(data = NA, ncol = 4)
> Dets <- vector()
> Traces <- vector()
> for (j in 1:Simulations) {
+     X <- rmnorm(nu, mean = c(0, 0), varcov = Sigma)
+     W <- t(X) %*% X
+     Dets <- c(Dets, det(W))
+     Traces <- c(Traces, sum(diag(W)))
+     W_xx <- c(W_xx, W[1, 1])
+     W_yy <- c(W_yy, W[2, 2])
+     W_xy <- c(W_xy, W[1, 2])
+     Vec_W <- rbind(Vec_W, matrix(W, nrow = 1, ncol = 4))
+ }
> Vec_W <- Vec_W[-1, ]

> library(matrixcalc)
> Expected_Value <- matrix(nu * Sigma, nrow = 1, ncol = 4)
> K_22 <- commutation.matrix(2, 2)
> Covariance <- nu * (diag(4) + K_22) %*% (Sigma %x% Sigma)
> Sample_Mean <- colMeans(Vec_W)
> Sample_Covariance <- cov(Vec_W)

> library(scatterplot3d)
> scatterplot3d(W_xx, W_xy, W_yy, highlight.3d = TRUE, col.axis = "blue",
+     col.grid = "lightblue", pch = 20, angle = 65)

> par(mfrow = c(2, 2))
> plot(W_xx, W_xy, pch = 19, col = "blue", main = "W_xy Vs W_xx")
> plot(W_yy, W_xy, pch = 19, col = "blue", main = "W_xy Vs W_yy")
> plot(W_xx, W_yy, pch = 19, col = "blue", main = "W_yy Vs W_xx")
> plot(W_xx * W_yy - W_xy * W_xy, pch = 19, col = "blue", main = "W_xx * W_xx - W_xy *W_xy")