Scienco.org

Although R supports a lot of different probability distributions, it does not support the scaled inverse chi-square distribution - well, not directly. But it can of course be implemented.

Let . Because of how the scaled inverse chi-square distribution is defined, it's possible to simulate from it in R as follows (checking the arguments is not included in the following code):

rinvchisq <- function(n, df, scale)
{
  return(scale * (1 / rchisq(n, df)))
}

When we're able to simulate from the distribution, we're able to find quantiles. But we can actually find those directly by using the definition of quantiles and the chi-squared distribution. Assume that we want to find the $\alpha$-quantile $x_\alpha$. Notice that with informal notation, the following holds:

If , then is the -quantile for . All together we now have:

Because of this, we can implement the following in R (checking the arguments is still not included in the following code):

qinvchisq <- function(p, df, scale)
{
  return(scale / qchisq(1 - p, df))
}

Checking the functions can be done like this:

invchisqcheck <- function(p, n, df, scale)
{
  cat("Check of rinvchisq and qinvchisq...\n")
  cat("Quantiles:", p, "\n")
  cat("Simulated with rinvchisq:\n")
  cat(quantile(rinvchisq(simulations, df, scale), quantiles), "\n")
  cat("Calculated with qinvchisq:\n")
  cat(qinvchisq(quantiles, df, scale), "\n")
}

invchisqcheck(c(0.025, 0.975), 100000, 10, 1)
invchisqcheck(c(0.025, 0.975), 100000, 11, 68)

For further information, please refer to http://en.wikipedia.org/wiki/Scale-inverse-chi-square_distribution.