3d-plot of Marchenko-Pastur law

marchenkopastur.R

# 3d-plot of Marchenko-Pastur law
# p 6 of High-Dim statistics (wainwright)
library(rgl)
alpha <- seq(0.01, 0.99, 0.01)
gamma <- seq(0, 4, length.out = 500)
fmp <- function(g, a) {
  gtmin <- g - (1 - sqrt(a))^2 
  tmaxg <- (1 + sqrt(a))^2 - g
  if (gtmin <= 0 | tmaxg <= 0) return(0)
  sqrt(tmaxg * gtmin) / g # see errata https://people.eecs.berkeley.edu/~wainwrig/highdim_errata.txt
}

fmpn <- function(g, a) {
  o <- Vectorize(fmp, "g")(g, a)
  o/sum(o)
}

persp3d(
  x = gamma, 
  y = alpha, 
  z = sapply(alpha, \(a) fmpn(gamma, a)), 
  theta = 110,
  expand = 0.5,
  zlab = "Density (rescaled)", 
  border = "#12888855", 
  col = "#12FAFA", 
  nticks = 10
)