## ----cranInstall, eval=FALSE--------------------------
# install.packages("invgamma")


## ----githubInstall, eval=FALSE------------------------
# remotes::install_github("dkahle/invgamma")


## ----load---------------------------------------------
library("invgamma")


## ----echo=FALSE, eval=FALSE---------------------------
# gamma(-3.5)
# pi/(sin(-3.5*pi)*gamma(4.5))


## ----echo=FALSE, eval=FALSE---------------------------
# gamma(3.2)
# 2.2*1.2*gamma(1.2)
# 
# gamma(-2.2)
# pi/(sin(-2.2*pi)*gamma(3.2))


## ----echo=FALSE, eval=FALSE---------------------------
# x <- 14.25
# lgamma(x)
# .5*log(2*pi) + (x-.5)*log(x) - x
# 
# gamma(x)
# sqrt(2*pi)*x^(x-.5)*exp(-x)


## -----------------------------------------------------
x <- 2.7; alpha <- 2; beta <- 3

dgamma(x, shape = alpha, scale = beta)

dpois(alpha-1, lambda = x/beta) / beta


## ----dinvgamma, fig.height=2.75, fig.cap="The PDF of the Inv-Gamma(7,10) distribution can be visualized with \\code{dinvgamma()}.", fig.pos="h!", message=FALSE, warning=FALSE----
# load helper packages, installing where necessary
install_if_needed <- function(pkg) {
  if (!requireNamespace(pkg)) install.packages(pkg)
}
install_if_needed("tidyverse"); library("tidyverse")
install_if_needed("patchwork"); library("patchwork")
library("scales", warn.conflicts = FALSE)
theme_set(theme_minimal())
theme_update(panel.grid.minor = element_blank())

# plot the inverse gamma PDF
ggplot() + 
  stat_function(
    fun = dinvgamma, n = 251, xlim = c(0, 6), alpha = .65,
    args = list(shape = 7, rate = 10), geom = "polygon"
  ) + theme(axis.title = element_blank())


## ----pinvgamma----------------------------------------
f <- function(x) dinvgamma(x, 7, 10)
q <- 2
integrate(f, 0, q)
(p <- pinvgamma(q, 7, 10))


## ----qinvgamma----------------------------------------
qinvgamma(p, 7, 10) # = q


## ----rinvgamma----------------------------------------
set.seed(1)
n <- 1e5
draws <- rinvgamma(n, 7, 10)
mean(draws <= q)


## ----correct, fig.height=3.25, fig.cap="(Left) A histogram based on 100,000 draws generated with \\code{rinvgamma()} superimposed with the inverse gamma density provided by \\code{dinvgamma()} (red); and (right) a QQ plot of the estimated quantiles (order statistics) vs theoretical quantiles computed with \\code{qinvgamma()}. Both strongly suggest the RNG is correct.", fig.pos="h!"----
df <- data.frame("x" = draws)

p_hist <- ggplot(df, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 250) +
  geom_function(fun = f,  color = "red", n = 251) + 
  coord_cartesian(xlim = c(0, 6))

qq_df <- data.frame(
  "theoretical_quantiles" = ppoints(n) |> qinvgamma(7, 10), 
  "observed_quantiles" = sort(draws)
)

p_qq <- ggplot(qq_df, aes(theoretical_quantiles, observed_quantiles)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  labs("x" = "Theoretical Quantiles", "y" = "Observed Quantiles")

p_hist + p_qq


## ----invgamma-rate-scale------------------------------
dinvgamma(.75, shape = 7, rate = 10)
dinvgamma(.75, shape = 7, scale = 1/10)


## ----warning=FALSE------------------------------------
test_rinvgamma <- function(shape, rate, n = 1e5) {
  sample <- rinvgamma(n, shape, rate)
  ks.test(sample, function(p) pinvgamma(p, shape, rate))$p.value
}

test_rinvgamma(shape = 3, rate = 7)


## ----gof-invgamma, fig.height=4, fig.cap="KS GoF tests of inverse gamma draws for different parameter values.", fig.pos="h!", warning=FALSE, message=FALSE, cache=TRUE----
# create the grid over the parameter space
n_grid <- 51
param_vals <- 10**seq(-4, 4, length.out = n_grid)
param_grid <- expand_grid("shape" = param_vals, "rate" = param_vals)

# run the simulations (in parallel)
install_if_needed("furrr"); library("furrr"); furrr_options(seed = TRUE)
install_if_needed("parallelly"); library("parallelly")
plan(multisession(workers = min(availableCores(), 10)))

param_grid <- param_grid |> 
  mutate("p_val" = future_map2_dbl(shape, rate, test_rinvgamma))


# plot the results
fmt <- scales::math_format(10**.x)
ggplot(param_grid, aes(shape, rate, color = p_val)) +
  geom_point() +
  scale_x_log10(expression(alpha), n.breaks = 10, labels = fmt(-5:5)) +
  scale_y_log10(expression(lambda), n.breaks = 10, labels = fmt(-5:5)) +
  scale_color_binned(breaks = c(0, .05, 1)) +
  labs(color = "p value") +
  coord_equal()


## ----a-warning----------------------------------------
rinvgamma(4, shape = .001, rate = 3)


## ----r-rgamma-problems--------------------------------
rgamma(10, shape = .001, rate = 3) == 0


## ----r-rgamma-problems-2------------------------------
rgamma(1e7, shape = 3, rate = 7) |> 
  table() |> sort(decreasing = TRUE) |> head(n = 5)


## ----gof-invchisq, fig.height=4, fig.cap="KS GoF tests of inverse chi-squared draws for different parameter values.", fig.pos="h!", warning=FALSE, message=FALSE, cache=TRUE----
test_rinvchisq <- function(df, ncp, n = 1e5) {
  sample <- rinvchisq(n, df, ncp)
  ks.test(sample, function(p) pinvchisq(p, df, ncp))$p.value
}

expand_grid("df" = param_vals, "ncp" = param_vals) |> 
  mutate("p_val" = future_map2_dbl(df, ncp, test_rinvchisq)) |> 
  ggplot(aes(df, ncp, color = p_val)) +
    geom_point() +
    scale_x_log10(expression(nu), n.breaks = 10, labels = fmt(-5:5)) +
    scale_y_log10("ncp", n.breaks = 10, labels = fmt(-5:5)) +
    scale_color_binned(breaks = c(0, .05, 1)) +
    labs(color = "p value") +
    coord_equal()


## ----gof-invexp, fig.height=1, fig.cap="KS GoF tests of inverse exponential draws for different parameter values.", fig.pos="h!", warning=FALSE, message=FALSE, cache=TRUE----
test_rinvexp <- function(rate, n = 1e5) {
  sample <- rinvexp(n, rate = rate)
  ks.test(sample, function(p) pinvexp(p, rate))$p.value
}

tibble("rate" = 10**seq(-4, 4, length.out = 2*n_grid)) |> 
  mutate("p_val" = future_map_dbl(rate, test_rinvexp)) |> 
  ggplot(aes(rate, 0, color = p_val)) +
    geom_point() +
    scale_x_log10(expression(lambda), n.breaks = 10, labels = fmt(-5:5)) +
    scale_color_binned(breaks = c(0, .05, 1), guide = FALSE) +
    theme(axis.text.y = element_blank(), axis.title.y = element_blank(),
          panel.grid.major.y = element_blank()) +
    coord_equal()