rinvgauss: The Shifted Inverse Gaussian Distribution

In rettopnivek/seqmodels: Distribution Functions for Sequential Sampling Models

Description

Random generation, density, distribution, and quantile functions for the shifted inverse gaussian (or Wald) distribution, parameterized for Brownian motion. kappa refers to the threshold, xi refers to the rate of evidence accumulation towards this threshold, tau is the shift in response times and sigma refers to the within-trial variability for the rate of evidence accumulation (the coefficient of drift, typically fixed to 1).

Usage

rinvgauss(n, kappa, xi, tau = as.numeric(c(0)), sigma = as.numeric(c(1)))

dinvgauss(
  t,
  kappa,
  xi,
  tau = as.numeric(c(0)),
  sigma = as.numeric(c(1)),
  ln = FALSE
)

pinvgauss(
  t,
  kappa,
  xi,
  tau = as.numeric(c(0)),
  sigma = as.numeric(c(1)),
  ln = FALSE,
  lower_tail = TRUE
)

qinvgauss(
  p,
  kappa,
  xi,
  tau = as.numeric(c(0)),
  sigma = as.numeric(c(1)),
  bounds = 3,
  em_stop = 20,
  err = 1e-08
)

minvgauss(kappa, xi, tau = as.numeric(c(0)), sigma = as.numeric(c(1)))

Arguments

n	the number of draws for random generation.
kappa	a vector of thresholds determining when a decision terminates (kappa > 0).
xi	a vector of drift rates, or rates of evidence accumulation (xi ≥ 0).
tau	a vector of shift parameters denoting the lowest possible time that can be observed (0 ≥ tau < t).
sigma	a vector of the within-trial variabilities (sigma > 0).
t	a vector of times ( t > 0 ).
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).
bounds	upper limit 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.

Value

dinvgauss gives the density, pinvgauss gives the distribution function, qinvgauss approximates the quantile function, minvgauss computes the descriptive moments (mean, variance, standard deviation, skew, and excess kurtosis), and rinvgauss generates random deviates.

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

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

Details

The inverse gaussian distribution describes the first passage times through a positive threshold kappa for a space and time homogenous Wiener diffusion process.

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.

References

Dagpunar, J. (1988). Principles of Random Variate Generation. Oxford: Clarendon Press.

Heathcote, A. (2004a). Fitting Wald and ex-Wald distributions to response time data: An example using functions for the S-PLUS package. Behavior Research Methods Instruments & Computers, 36, 678 - 694.

Heathcote, A. (2004b). rtfit.ssc. Retrieved May 5, 2017 from Psychonomic Society Web Archive: http://www.psychonomic.org/ARCHIVE/.

Examples

# Density
dinvgauss( .9758, kappa = 1.0, xi = 1.0, tau = 0.3 )
# Distribution function
pinvgauss( .9758, kappa = 1.0, xi = 1.0, tau = 0.3 )
# Quantile function (Accurate to ~4 decimal places)
round( qinvgauss( p = .5, kappa = 1.0, xi = 1.0, tau = 0.3 ), 4 )
# Descriptive moments
minvgauss( kappa = 1.0, xi = 1.0, tau = 0.3 )

# Simulation (No shift)
sim <- rinvgauss( 1000, kappa = 0.8, xi = 2.0 )

# Function to obtain maximum likelihood estimates
param_est <- function( dat, tau_hat = 0 ) {
  # Estimate threshold and drift from first two moments of
  # data (Heathcote, 2004):
  dat <- dat - tau_hat # Apply shift
  xi_hat <- sqrt( mean( dat )/var( dat ) );
  kappa_hat <- xi_hat * mean( dat );
  return( c( kappa = kappa_hat, xi = xi_hat, tau = tau_hat ) )
}
print( param_est( sim ) )

# Non-zero shift parameter
sim <- rinvgauss( 1000, kappa = 1.6, xi = 1.5, tau = .3 )

# Estimating shift parameter

# Function to compute sum of log-likelihoods
f <- function( tau_hat, dat ) {
  prm = param_est( dat, tau_hat = tau_hat )
  sll = sum( dinvgauss( dat, prm[1], prm[2], tau = prm[3], ln = T ) )
  return( sll )
}
tau_hat <- optimize( f, c( 0.0, min( sim ) ), dat = sim, maximum = T )
print( param_est( sim, tau_hat = tau_hat$maximum ) )

# Plotting
layout( matrix( 1:4, 2, 2, byrow = T ) )
# Parameters
prm <- c( k = 0.8, x = 1.6, t = 0.3, s = 1.0 )
# Density
obj <- quickdist( 'sig', 'PDF', prm )
plot( obj ); lines( obj )
# CDF
obj <- quickdist( 'sig', 'CDF', prm )
plot( obj ); lines( obj )
# Quantiles
obj <- quickdist( 'sig', '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( 'sig', 'HF', prm )
plot( obj ); lines( obj )

rettopnivek/seqmodels documentation built on May 1, 2020, 2:59 p.m.