inverse gamma (0.001,0.001) prior on the variance in the Bayesian hierarchical model This 8 schools data is from Gelman 2006 paper: http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf. In Figure 1 (c), the prior density of inverse gamma (0.001,0.001) was overlain on the posterior histogram. The author claims that "This prior distribution is even more sharply peaked near zero and further distorts posterior inferences", which I find difficult to understand because inverse gamma (0.001,0.001) only has a very small probability near 0, e.g., $P(\sigma < 1 | \alpha =\beta=.001) = 0.006$ for the prior but the posterior $P(\sigma < 1) = .487$. So how does this prior distort the posterior towards 0?

UPDATE: Thanks so much for your reproducible code. I have accepted it as an answer.  It appears that the prior density in the plot has been scaled to the range (0,30), that's why we see a spike of inv-Gamma(0.001, 0.001) near 0. But the issue is that if you draw samples from inv-Gamma(0.001, 0.001), you are going to see the probability $\sigma < 1$ is quite small and most values would be very large, so I do not understand how this prior can still distort the posterior to have a spike near 0. PS: the CDF of the inverse gamma distribution on Wiki page is defined using UPPER incomplete gamma function, so you should look at the "uppinc" when you run pracma::gammainc(shape, scale / 1) / gamma(shape)
 A: The point here is that the posterior for a invgamma(1, 1) prior has greater density near zero than that for a Uniform prior, and invgamma(.001, .001) has even more density near zero than invgamma(1, 1). In other words, you should compare the posterior distribution for each plot, and compare the prior distributions for each plot, but not necessarily compare the prior to the posterior within each plot.
A: You're right that the $\text{inverse-gamma}(0.001, 0.001)$ prior assigns a small probability to values of $\sigma$ near $0$. As you noted, $P(\sigma < 1 | \alpha =\beta=.001) = 0.006$. However, relative to other values of $\sigma$, the most likely value of $\sigma$ is near $0$. The density near 0 is high relative to other regions. While the density assigns low absolute probability to $\sigma$ near 0, the most likely value is $\sigma$ is still near 0.
Contrast this to the $\text{inverse-gamma}(1, 1)$ prior. This prior still says that $\sigma$ should be close to 0. However, the tails of this prior drop more slowly. This means larger values of $\sigma$ are still considered plausible. Contrast both inverse-gamma priors to the $\text{Unif}[0, A]$ prior which gives equal likelihood to small and large values of $\sigma$.
The author is arguing that the inverse-gamma distribution is a bad choice for a noninformative prior. It puts most of it's density near 0, and therefore pulls the posterior distribution close to 0. Using this class of priors is like saying you expect a small standard deviation before you see any data. This is not the behavior people usually associate with non-informative priors. You're right that the posterior isn't as peaked near 0 as the prior. However, it would be even further from 0 if the prior didn't assign so much mass at 0. The $Unif[0, A]$ prior doesn't pull the posterior to 0, and, instead, let's it spread out further of the real number line.
A: Figure 1 in the "Priors for variances" paper compares three prior distributions for the hierarchical standard deviation, $\sigma_\alpha$, in a two-level normal hierarchical model. It illustrates the drawbacks of the supposedly non-informative inverse gamma prior on the eight schools example.
As the answers by @Eoin and @Eli point out, the three priors are ordered so that, from left to right, each prior puts more probability mass near zero. In practice the inverse gamma priors suggest that we expect the school effects $\alpha_j$ shrink towards 0 even before looking at the data. In other words, the uniform prior is the least informative and the inverse-gamma($\epsilon$, $\epsilon$) prior is the most informative.
The informativeness of priors is to some extent defined with respect to the likelihood: a prior is strongly informative if it doesn't support values that are consistent with the data. For the eight schools, the likelihood is close to flat in the range [0,10], so any prior that concentrates only on a part of this range ends up being much more informative than the likelihood and the data.
This is perhaps easier to see in the plots below which show the prior (red curve), the likelihood (blue curve) and the posterior (histogram) because
$$
posterior \propto prior \times likelihood
$$
With the uniform prior, the posterior is informed mostly by the likelihood, ie, by the data. But with the inverse gamma prior, the posterior doesn't look like much different from the prior even after observing the data from the eight schools. Why bother with collecting data if we are not prepared to update our prior in the face of evidence?
Note: The likelihood for eight schools is a function of both the overall effect $\mu$ and the hierarchical standard deviation $\sigma_\alpha$. I picked $\mu$ = 10 to plot the likelihood as a function of $\sigma_\alpha$ but the choice doesn't matter: for any value of $\mu$ in [0,15] the likelihood supports values for $\sigma_\alpha$ in the range [0,10].

Finally, the inverse-gamma($\epsilon$, $\epsilon$) has most of its probability mass close to zero. It seems you have inadvertently computed the upper-tail probability P{X > x} instead of the lower-tail probability P{X ≤ x}.
I'm sure you've looked up the definition: The CDF of the inverse gamma distribution with parameters shape $\alpha$ and scale $\beta$ is $F(x;\alpha,\beta) = \Gamma(\alpha,\beta/x) / \Gamma(\alpha)$. The numerator is the incomplete gamma function, the numerator is the (complete) gamma function.
Let's do the computation in R; we will use the pracma::gammainc function which helpfully returns both the upper and lower incomplete gammas, so any confusion is avoided.
shape <- scale <- 0.001

# lowinc = P{X ≤ x}; uppinc = P{X > x} for x = 0.1 and x = 1

pracma::gammainc(shape, scale / 0.1) / gamma(shape)
#>       lowinc       uppinc       reginc 
#> 0.0933783135 0.0061116005 0.0009391118
pracma::gammainc(shape, scale / 1) / gamma(shape)
#>       lowinc       uppinc       reginc 
#> 0.9936876467 0.0063123533 0.0009942606


The code to reproduce the figures, in all its gory details:
library("tidyverse")
library("brms")

set.seed(2022)

y_j <- c(28, 8, -3, 7, -1, 1, 18, 12)
sigma_j <- c(15, 10, 16, 11, 9, 11, 10, 18)

schools <- tibble(
  school = letters[1:8], y_j, sigma_j
)

log_likelihood <- function(mu, sigma_a) {
  sum(dnorm(y_j, mean = mu, sd = sqrt(sigma_j^2 + sigma_a^2), log = TRUE))
}
likelihood <- function(mu, sigma_a) {
  exp(log_likelihood(mu, sigma_a))
}

scaled_curve <- function(func, from = 0, to = 1, n = 1001, ...) {
  x <- seq(from = from, to = to, length.out = n)
  y <- eval(substitute(func), envir = list(x = x))
  # Scale to match a density histogram
  y <- y / sum(y, na.rm = TRUE) / ((to - from) / n) * .95
  lines(x, y, ...)
}

schools_posterior <- function(sigma_a_prior) {
  fit <- brm(
    y_j | se(sigma_j) ~ 1 + (1 | school),
    family = gaussian,
    prior = prior_string(sigma_a_prior, class = "sd"),
    data = schools
  )
  hist(
    as_draws_array(fit, variable = "sd_school__Intercept"),
    breaks = 60, probability = TRUE, col = NULL,
    xlim = c(0, 30), xlab = expression(sigma[alpha]),
    yaxt = "n", ylab = NULL,
    main = paste("The prior is", sigma_a_prior)
  )
  scaled_curve(
    sapply(x, function(x) likelihood(10, x)),
    from = 0, to = 30, col = "blue", lwd = 2
  )
}

schools_posterior(
  "uniform(0,100)"
)
scaled_curve(
  dunif(x, min = 0, max = 100),
  from = 0, to = 30, col = "red", lwd = 2
)

schools_posterior(
  "inv_gamma(1,1)"
)
scaled_curve(
  MCMCpack::dinvgamma(x, 1, 1),
  from = 0, to = 30, col = "red", lwd = 2
)

schools_posterior(
  "inv_gamma(0.001,0.001)"
)
scaled_curve(
  MCMCpack::dinvgamma(x, 0.001, 0.001),
  from = 0, to = 30, col = "red", lwd = 2
)

