
Sampling from Gamma Distribution
srgamma_custom.Rd
The srgamma_custom()
function generates random samples from a Gamma distribution using the STORS algorithm.
It employs an optimized proposal distribution around the mode and Adaptive Rejection Sampling (ARS) for the tails.
Value
A numeric vector of length n
containing random samples from the Gamma distribution.
The shape
and rate
parameters are specified during the optimization process using srgamma_optimize()
.
NOTE: When the x
parameter is specified, it is updated in-place with the simulation for performance reasons.
Details
The Gamma Distribution
The Gamma distribution has the probability density function (PDF): $$f(x | \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} \exp(-\beta x), \quad x \geq 0,$$ where:
- \(\alpha\)
is the shape parameter (\(\alpha > 0\)), which determines the shape of the distribution.
- \(\beta\)
is the rate parameter (\(\beta > 0\)), which determines the rate of decay.
The Gamma distribution is widely used in statistics, particularly in Bayesian inference and modelling waiting times.
This function samples from a proposal constructed using srgamma_optimize
, employing the STORS algorithm.
By default, srgamma_custom()
samples from the standard Gamma distribution with shape = 1
and rate = 1
.
The proposal distribution is pre-optimized at package load time using srgamma_optimize()
with
steps = 4091
, creating a scalable proposal centred around the mode.
Note
This function is not scalable. Therefore, only the srgamma_custom()
version is available, which requires the proposal to be pre-optimized using srgamma_optimize()
before calling this function.
See also
srgamma_optimize
to optimize the custom proposal.
Examples
# Generate 10 samples from Gamma Distribution
samples <- srgamma_custom(10)
print(samples)
#> [1] 0.952679445 1.254666123 1.144112999 0.001567012 0.251546634 0.963247523
#> [7] 1.729025469 1.781825998 1.186380538 0.407005608
# Generate 10 samples using a pre-allocated vector
x <- numeric(10)
srgamma_custom(10, x = x)
#> [1] 2.4761374 0.3352764 0.1282204 1.3389049 0.2770518 2.5073452 0.2811029
#> [8] 2.0085429 0.1123188 1.1362554
print(x)
#> [1] 2.4761374 0.3352764 0.1282204 1.3389049 0.2770518 2.5073452 0.2811029
#> [8] 2.0085429 0.1123188 1.1362554