Findalphabeta_gamma: Quantile difference when tuning Gamma distribution...
In andrewzm/bicon: Bivariate modelling using the conditional approach

Description Usage Arguments Details Value Examples

See details below.

1	Findalphabeta_gamma(pars, p5, p95)

`pars`	the shape and rate parameters (in that order) of the Gamma distribution
`p5`	the 5-th percentile of the desired error distribution
`p95`	the 95-th percentile of the desired error distribution

Assume that in our model there is an error term (for example, observation error) e \sim \mathcal{N}(0,σ^2) and we are required to estimate the precision, σ^{-2}. Assume further that we have sufficient prior knowledge to know that F_e(0.05) = p5 and F_e(0.95) = p95, where F_e is the cumulative distribution function of e and p5 and p95 are known. This function helps configure a prior distribution over σ^{-2} which reflects this prior belief. Specifically, let Λ_{α,β} be the cumulative distribution function of the Gamma distribution with parameters α and β. Then this function returns the discrepancy

\frac{Λ_{α,β}(0.05) - 1/(p95)^2}{1/(p95)^2} + \frac{Λ_{α,β}(0.95) - 1/(p5)^2}{1/(p5)^2}

For example, we can see how representative a Gamma distribution over the precision with shape and rate parameters in pars is of our prior belief that p5 and p95 are 1 and 5, respectively, by calling Findalphabeta_gamma(pars,1,5). This function can hance be passed to an optimisation routine to find better values for pars.

A measure of discrepancy between the Gamma distribution with parameters par and the desired 5-th 95-th percentiles of the error signal.

require(actuar)
# Find the Gamma distribution over the precision corresponding to the
# prior belief of the error (1/sqrt(precision)) lying between p5=2, p95=5
initpars <- c(5,0.1)
hyp_pars <- optim(par=initpars,Findalphabeta_gamma, p5=1, p95=5)

# Now simulate from a Gamma with these parameters and verify quantiles
X <- rgamma(shape = hyp_pars$par[1], rate = hyp_pars$par[2],n=10000)
print( quantile(1/sqrt(X),c(0.05,0.95)))