remg: The Exponentially Modified Gaussian Distribution
In rettopnivek/seqmodels: Distribution Functions for Sequential Sampling Models

Description Usage Arguments Value Details References Examples

View source: R/RcppExports.R

Random generation, density, distribution, quantile, and descriptive moments functions for the the convolution of gaussian and exponential random variables (the ex-gaussian distribution), with location equal to mu, scale equal to sigma and rate equal to lambda.

remg(n, mu, sigma, lambda)

demg(x, mu, sigma, lambda, ln = FALSE)

pemg(q, mu, sigma, lambda, ln = FALSE, lower_tail = TRUE)

qemg(
  p,
  mu,
  sigma,
  lambda,
  bounds = as.numeric(c(0, 3)),
  em_stop = 20,
  err = 1e-08
)

memg(mu, sigma, lambda)

`n`	number of observations to be generated.
`mu`	vector of location parameters for the gaussian variable.
`sigma`	vector of scale parameters (i.e., standard deviations) for the gaussian variable (sigma > 0).
`lambda`	vector of rate parameters for the exponential variable (lambda > 0).
`x, q`	vector of quantiles.
`ln`	logical; if `TRUE`, probabilities are given as log(p).
`lower_tail`	logical; if `TRUE` (default), probabilities are P(X ≤ x) otherwise P( X > x).
`p`	vector of probabilities.
`bounds`	lower and upper limits of the quantiles to explore for the approximation via linear interpolation.
`em_stop`	the maximum number of iterations to attempt to find the quantile via linear interpolation.
`err`	the number of decimals places to approximate the cumulative probability during estimation of the quantile function.

demg gives the density, pemg gives the distribution function, qemg gives the quantile function, memg computes the descriptive moments (mean, variance, standard deviation, skew, and excess kurtosis), and remg generates random deviates.

The length of the result is determined by n for remg, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result.

An exponentially modified gaussian distribution describes the sum of independent normal and exponential random variables, possessing a characteristic positive skew due to the exponential variable. The ex-gaussian distribution is therefore useful for providing descriptive fits to response time distributions (e.g., Luce, 1986).

A linear interpolation approach is used to approximate the quantile function, estimating the inverse of the cumulative distribution function via an iterative procedure. When the precision of this estimate is set to 8 decimal places, the approximation will be typically accurate to about half of a millisecond.

The example section demonstrates how to compute maximum likelihood estimates based on the moments from a set of data.

Luce, R. D. (1986). Response times: Their role in inferring elementary mental organization. New York, New York: Oxford University Press.

Terriberry, T. B. (2007). Computing Higher-Order Moments Online. Retrieved from https://people.xiph.org/~tterribe/notes/homs.html

# Density function
demg( x = 0.8758, mu = 0.0, sigma = 1.0, lambda = 1.0 )
# Distribution function
pemg( q = 0.8758, mu = 0.0, sigma = 1.0, lambda = 1.0 )
# Quantile function (Accurate to ~4 decimal places)
round( qemg( p = .5, mu = 0.0, sigma = 1.0, lambda = 1.0 ), 4 )
# Descriptive moments
memg( mu = 0.44, sigma = 0.07, lambda = 2.0 )

# Simulation
sim <- remg( n = 1000, mu = 0.44, sigma = 0.07, lambda = 2.0 );

# Function to obtain maximum likelihood estimates
param_est <- function( dat ) {
  # Compute 1st, 2nd, and 3rd moments ( Terriberry, 2007).
  n <- 0; m <- 0; m2 <- 0; m3 <- 0;
  for ( k in 1:length( dat ) ) {
    n1 <- n; n <- n + 1
    term1 <- dat[k] - m; term2 = term1/ n; term3 = term1 * term2 * n1
    # Update first moment
    m <- m + term2
    # Update third moment
    m3 <- m3 + term3 + term2 * (n - 2) - 3 * term2 * m2
    # Update second moment
    m2 <- m2 + term3
  }
  # Compute standard deviation of sample
  s <- sqrt( m2 / ( n - 1.0 ) )
  # Compute skewness of sample
  y <- sqrt( n ) * m3 / ( m2^1.5 );
  # Estimate parameters
  mu_hat <- m - s * ( y/2 )^(1/3)
  sigma_hat <- sqrt( (s^2) * ( 1 - (y/2)^(2/3) ) )
  lambda_hat <- 1/( s * (y/2)^(1/3) )
  return( c( mu = mu_hat, sigma = sigma_hat, lambda = lambda_hat ) )
}

print( param_est( sim ) )

# Plotting
layout( matrix( 1:4, 2, 2, byrow = T ) )
# Parameters
prm <- c( m = .44, s = .07, l = 2.0 )
# Density
obj <- quickdist( 'emg', 'PDF', prm )
plot( obj ); lines( obj )
# CDF
obj <- quickdist( 'emg', 'CDF', prm )
plot( obj ); lines( obj )
# Quantiles
obj <- quickdist( 'emg', 'QF', prm, x = seq( .2, .8, .2 ) )
plot( obj ); prb = seq( .2, .8, .2 );
abline( h = prb, lty = 2 );
lines( obj, type = 'b', pch = 19 )
# Hazard function
obj <- quickdist( 'emg', 'HF', prm )
plot( obj ); lines( obj )