sample.distr: Randomly sample from a user defined probability distribution

In Chitran1987/StatsChitran: StatsChitran

Randomly sample from a user defined probability distribution

Description

Returns a set of values of the random variable, \bar{X}, sampled from a user defined probability distribution f(\bar{X}, \hat{\theta}).
Identical functionality to the runif and rnorm functions present in the base R package.

Usage

sample.distr(n, func, f.par, xmin, xmax)

Arguments

n	The no. of values of \bar{X} that are expected to be returned. Will be the length of the return vector.
func	The function f(\bar{X}, \hat{\theta}). The argument passed to func must be a function. The first argument to the function that is passed to func must be the unknown \bar{X} and NOT the distribution parameters, \hat{\theta}. The subsequent arguments must encode the distribution parameters \hat{\theta} as the user desires. The ordering of the arguments for the function passed to func is important. Please see examples.
f.par	The specific values of the parameters, \hat{\theta}, that define the distribution, f(\bar{X}, \hat{\theta}), in this specific case. Must be a numeric atomic vector. It is important that the length of f.par is one less than the no. of arguments of the function passed to func. This is because the first argument to func would be \bar{X}. Ordering of the values of the elements of f.par should adhere to the ordering of the parameters, \hat{\theta}, in f(\bar{X}, \hat{\theta}). Please see examples.
xmin, xmax	Numeric values within which the returned values of \bar{X} will lie

Details

Uses rejection sampling to arrive at the results of \bar{X}.
Thus results for n > 10^7 could strain system resources, since rejection sampling is a monte-carlo method and needs overgeneration beyond n in the rejection space.
The range of the function f(\bar{X},\hat{\theta}) input to func will be normalized to f(\bar{X}, \hat{\theta}) \in [0, 1].
Hence it should be ensured that f(\bar{X},\hat{\theta}) is devoid of badly behaved singularities within \code{xmin} \leq \bar{X} \leq \code{xmax}.

Value

The reurned value is a 2 element list

• random_var: The random variable \bar{X} which follows the distribution f(\bar{X}, \hat{\theta}) returned as a R vector

• likelihood: The density value f(\bar{X}, \hat{\theta}) for each corresponding value in \bar{X} as a R vector

Author(s)

Chitran Ghosal

References

https://en.wikipedia.org/wiki/Rejection_sampling

Examples

####Define the probability distribution
sin.exp <- function(X, al, omeg){
  y <- exp(-al*abs(X))*(sin(omeg*X))^2
  return(y)
}

####define the distribution parameters
par_f <- c(2, 2)

####plot the bare distribution function
RV <- seq(-10, 10, by=0.01) #Random variable needed for the plotting
Y <- nrm(sin.exp(X=RV, al =par_f[1], omeg = par_f[2])) #The density/likelihood values against the given random variables
plot(RV, Y, type = 'l')

####retrieve random values of X out of the distribution and plot its histogram for different values of n

## for n = 10
L <- sample.distr(n = 10^1, func = sin.exp, f.par = par_f, xmin = -10, xmax = 10)
h <- hist(L$random_var, probability = T, breaks = seq(-10, 10, by = 0.25), border = NA, main = "### n = 10 ")
lines(RV, nrm(Y, min = 0, max = max(h$density)), col = 'black', lwd=2.5)

## for n = 100
L <- sample.distr(n = 10^2, func = sin.exp, f.par = par_f, xmin = -10, xmax = 10)
h <- hist(L$random_var, probability = T, breaks = seq(-10, 10, by = 0.25), border = NA, main = "### n = 100 ")
lines(RV, nrm(Y, min = 0, max = max(h$density)), col = 'black', lwd=2.5)

## for n = 1000
L <- sample.distr(n = 10^3, func = sin.exp, f.par = par_f, xmin = -10, xmax = 10)
h <- hist(L$random_var, probability = T, breaks = seq(-10, 10, by = 0.25), border = NA, main = "### n = 1000 ")
lines(RV, nrm(Y, min = 0, max = max(h$density)), col = 'black', lwd=2.5)

## for n = 10^4
L <- sample.distr(n = 10^4, func = sin.exp, f.par = par_f, xmin = -10, xmax = 10)
h <- hist(L$random_var, probability = T, breaks = seq(-10, 10, by = 0.25), border = NA, main = "### n = 10000 ")
lines(RV, nrm(Y, min = 0, max = max(h$density)), col = 'black', lwd=2.5)

Chitran1987/StatsChitran documentation built on Feb. 23, 2025, 8:30 p.m.