# Ox Standard Packages Reference

## Chapter contents:

| Arma |
arma0 | armaforc | armagen | armavar
diffpow | modelforc | pacf
| Maximization |
GetMaxControl | GetMaxControlEps | MaxBFGS | MaxControl
MaxControlEps | MaxConvergenceMsg | MaxFSQP | MaxNewton
MaxSimplex | Num1Derivative | Num2Derivative | NumJacobian
SolveQP
| Probability |
densbeta | densgamma | densgh | densgig | denskernel
densmises | denspoisson
probbeta | probbvn | probgamma | probmises | probmvn
probpoisson | quanbeta | quangamma | quanmises
ranbeta | ranbinomial | ranbrownianmotion | ranchi | randirichlet
ranexp | ranf | rangamma | rangh | rangig
ranlogarithmic | ranlogistic | ranlogn | ranmises | ranmultinomial
rannegbin | ranpoisson | ranpoissonprocess | ranshuffle | ranstable
ransubsample | rant | ranuorder | ranwishart
| QuadPack |
QAG | QNG | QAGS | QAGP | QAGI
QPEPS

## Arma package

### arma0

    #include <arma.h>
    arma0(const ma, const vp, const cp, const cq);

`ma`
  in: T `x` n matrix A
`vp`
  in: 1 `x` s matrix with AR coefficients \phi1, \phi2, ..., \phip followed by MA coefficients \theta1, \theta2, ..., \thetaq, s >= p+q
`cp`
  in: int, no of AR coefficients (could be 0)
`cq`
  in: int, no of MA coefficients (could be 0)

Return value
  Returns the residual from applying the ARMA(p,q) filter to each column of A. The result has the same dimensions as `ma`. The first p rows of the return value will be zero.

  For example when p = 1 and q = 2:

    e0 = 0,
    e1 = a1 - \phi1 a0 - \theta1 e0,
    e2 = a2 - \phi1 a1 - \theta1 e1 - \theta 1 e0, etc.
See also
  `armagen`, `armaforc`, `armavar`, `diff0`, `diffpow`, `pacf`

### armaforc

    #include <arma.h>
    armaforc(const mx, const vp, const cp, const cq);
    armaforc(const mx, const vp, const cp, const cq,
        const ma);
    armaforc(const mx, const vp, const cp, const cq,
        const ma, const me);

`mx`
  in: H `x` n matrix X, fixed part of forecasts
`vp`
  in: 1 `x` s matrix with AR coefficients \phi1, \phi2, ..., \phip followed by MA coefficients \theta1, \theta2, ..., \thetaq, s >= p+q
`cp`
  in: int, no of AR coefficients (could be 0)
`cq`
  in: int, no of MA coefficients (could be 0)
`ma`
  in: (optional argument) T `x` n matrix A, pre-forecast data values (default is zero)
`me`
  in: (optional argument) T `x` n matrix E, pre-forecast residual values (default is zero)

Return value
  Returns the forecasts from an ARMA(p,q) model as an H `x` n matrix. The same model is applied to each column of `mx`.

  For example when p = 2 and q = 2:

    ahat0 = x0 + \phi1 aT-1 + \phi2 aT-2 + \theta1 eT-1 + \theta2 eT-2,
    ahat1 = x1 + \phi1 ahat0 + \phi2 aT-1 + \theta2 eT-1,
    ahat2 = x2 + \phi1 ahat1 + \phi2 ahat0, etc.
See also
  `arma0`, `armavar`, `cumulate`, `modelforc`

### armagen

    #include <arma.h>
    armagen(const mx, const me, const vp, const cp, const cq);

`mx`
  in: T `x` n matrix of known components X
`me`
  in: T `x` n matrix of errors E
`vp`
  in: 1 `x` s matrix with AR coefficients \phi1, \phi2, ..., \phip followed by MA coefficients \theta1, \theta2, ..., \thetaq, s >= p+q
`cp`
  in: int, no of AR coefficients (could be 0)
`cq`
  in: int, no of MA coefficients (could be 0)

Return value
  Generates a an ARMA(p,q) series from an error term (`me`) and a mean term (`mx`) The result has the same dimensions as `mx`. The first p rows of the return value will be identical to those of `mx`; the recursion will be applied from the pth term onward (missing lagged errors are set to zero).

  For example when p = 1 and q = 2:

    a0 = x0,
    a1 = x1 + \phi1 a0 + e1 + \theta1 e0,
    a2 = x2 + \phi1 a1 + e2 + \theta1 e1 + \theta2 e0, etc.
See also
  `arma0`, `armaforc`, `armavar`, `cumsum`, `cumulate`

### armavar

    #include <arma.h>
    armavar(const vp, const cp, const cq, const dvar,
        const ct);

`vp`
  in: 1 `x` s matrix with AR coefficients \phi1, \phi2, ..., \phip followed by MA coefficients \theta1, \theta2, ..., \thetaq, s >= p+q
`cp`
  in: int, no of AR coefficients (could be 0)
`cq`
  in: int, no of MA coefficients (could be 0)
`dvar`
  in: double, variance of disturbance, \sigma2 \epsilon.
`ct`
  in: int, number of autocovariance terms required, T

Return value
  Returns a 1 `x` T matrix with the autocovariances of the ARMA(p,q) process. Or `0` if the computations failed (e.g. when all AR coefficients are zero).
See also
  `arma0`, `pacf`

### diffpow

    #include <arma.h>
    diffpow(const ma, const d);

`ma`
  in: T `x` n matrix A
`d`
  in: double, length of difference d, |d| <= 10000

Return value
  Returns a T `x` n matrix with (1-L)dA. The result has the same dimensions as `ma`.
See also
  `arma0`, `diff0`

### modelforc

    #include <arma.h>
    modelforc(const mU, const mData, const miDep,
        const miSel, const miLag, const mPi, const iTmin);

`mU`
  in: 0 or (T-T0) `x` n matrix with error terms
`mData`
  in: T `x` d matrix, the database
`miDep`
  in: 1 `x` n matrix with database indices of dependent variables
`miSel`
  in: 1 `x` k matrix with database indices of explanatory variables
`miLag`
  in: 1 `x` k matrix with lag lengths of explanatory variables
`miPi`
  in: n `x` k matrix with coefficients
`iTmin`
  in: int, first forecast observation (may be 0)

Return value
  Returns the dynamic forecasts from a linear dynamic model as a (T-T0) `x` n matrix.
See also
  `armaforc` `cumulate`

### pacf

    #include <arma.h>
    pacf(const macf);
    pacf(const macf, const alogdet);
    pacf(const macf, const alogdet, const my);
    pacf(const macf, const meps);

`macf`
  in: arithmetic type, T `x` 1 matrix of autocovariances or autocorrelations
`alogdet`
  in: (optional argument) address of variable
  out: the determinant of the filter
`my`
  in: (optional argument) T `x` n data matrix to apply filter to
`meps`
  in: (optional argument) T `x` n data matrix to apply inverse filter to

Return value
  - `pacf(macf);`
  - `pacf(macf, alogdet);`
    Returns a T `x` 1 matrix with the partial autocorrelation function of the columns of `macf`.

  - `pacf(macf, alogdet, my);`
    Returns a T `x` (n+1) matrix with the residuals from the filter based on the specified ACF applied to the columns of `my`. The last column contains the standard devations of the filter.

  - `pacf(macf, meps);`
    Returns a T `x` n matrix with the generated data from the inverse filter based on the specified ACF applied to the columns of `meps`.
  Returns 0 if the computations fail (the stochastic process has a root on the unit circle).
See also
  `acf`, `arma0`, `armavar`, `solvetoeplitz`

## Maximization package

### GetMaxControl, GetMaxControlEps

    #import <maximize>
    GetMaxControl();
    GetMaxControlEps();

Return value
  Return an array with three values and two values respectively.
  `GetMaxControl` returns `{ mxIter, iPrint, bCompact }`.
  `GetMaxControlEps` returns `{ dEps1, dEps2 }`.
See also
  `MaxControl`, `MaxControlEps`

### MaxBFGS

    #import <maximize>
    MaxBFGS(const func, const avP, const adFunc,
        const amInvHess, const fNumDer);

`func`
  in: a function computing the function value, optionally with derivatives
`avP`
  in: address of p `x` 1 matrix with starting values
  out: p `x` 1 matrix with final coefficients
`adFunc`
  in: address
  out: double, final function value
`amInvHess`
  in: `0` or address of p `x` p matrix, initial (inverse negative) quasi-Hessian K; a possible starting value is the identity matrix (which is used when `0` is used as argument);
  out: if not `0` on input: final K (not reliable as estimate of actual Hessian)
`fNumDer`
  in: 0: `func` provides analytical first derivatives 1: use numerical first derivatives

The supplied `func` argument should have the following format:

    func(const vP, const adFunc, const avScore, const amHessian);

`vP`
  in: p `x` 1 matrix with coefficients
`adFunc`
  in: address
  out: double, function value at `vP`
`avScore`
  in: 0, or an address
  out: if !0 on input: p `x` 1 matrix with first derivatives at `vP`
`amHessian`
  in: always 0 for MaxBFGS, as it does not need the Hessian;

returns
  1: successful, 0: function evaluation failed

Return value
  Returns the status of the iterative process:
  `MAX_CONV`
    Strong convergence Both convergence tests were passed, using tolerance eps = eps1.
  `MAX_WEAK_CONV`
    Weak convergence (no improvement in line search) The step length si has become too small. The convergence test ) was passed, using tolerance eps = eps2.
  `MAX_MAXIT`
    No convergence (maximum no of iterations reached)
  `MAX_LINE_FAIL`
    No convergence (no improvement in line search) The step length si has become too small. The convergence test was not passed, using tolerance eps = eps2.
  `MAX_FUNC_FAIL`
    No convergence (function evaluation failed)
  `MAX_NOCONV`
    No convergence Probably not yet attempted to maximize.

  The chosen default values for the tolerances are:

  eps1=1e-4, eps2=5e-3.

See also
  `MaxControl`, `MaxConvergenceMsg`, `Num1Derivative`, `Num2Derivative`

### MaxControl, MaxControlEps

    #import <maximize>
    MaxControl(const mxIter, const iPrint);
    MaxControl(const mxIter, const iPrint, const bCompact);
    MaxControlEps(const dEps1, const dEps2);

`mxIter`
  in: int, maximum number of iterations; default is 1000, use -1 to leave the current value unchanged
`iPrint`
  in: int, print results every iPrint'th iteration; default is 0, use -1 to leave the current value unchanged
`bCompact`
  in: int, if TRUE uses compact format for iteration results (optional argument)
`dEps1`
  in: double, eps1, default is 1e-4, use 0 to leave the current value unchanged
`dEps2`
  in: double, eps2, default is 5e-3, use 0 to leave the current value unchanged

Return value
  No return value.
See also
  `GetMaxControl`, `GetMaxControlEps`, `MaxBFGS`, `MaxNewton`, `MaxSimplex`

### MaxConvergenceMsg

    #import <maximize>
    MaxConvergenceMsg(const iCode);

`iCode`
  in: int, code returned by `MaxBFGS` or `MaxSimplex`

Return value
  Returns the text corresponding to the convergence code listed under the return values of `MaxBFGS`.
See also
  `MaxBFGS`, `MaxSimplex`

### MaxFSQP

`MaxFSQP` implements a feasible sequential quadratic programming technique to maximize a non-linear function subject to non-linear constraints. Unfortunately, this function was not ready in time for this release.

### MaxNewton

    #import <maximize>
    MaxNewton(const func, const avP, const adFunc,
        const amHessian, const fNumDer);

`func`
  in: a function computing the function value, optionally with derivatives
`avP`
  in: address of p `x` 1 matrix with starting values
  out: p `x` 1 matrix with final coefficients
`adFunc`
  in: address
  out: double, final function value
`amHessian`
  in: `0` or address
  out: if not `0` on input: final Hessian matrix H
`fNumDer`
  in: 0: `func` provides analytical second derivatives 1: use numerical second derivatives

The supplied `func` argument should have the following format:

    func(const vP, const adFunc, const avScore, const amHessian);

`vP`
  in: p `x` 1 matrix with coefficients
`adFunc`
  in: address
  out: double, function value at `vP`
`avScore`
  in: 0, or an address
  out: if !0 on input: p `x` 1 matrix with first derivatives at `vP`
`amHessian`
  in: 0, or an address
  out: if !0 on input: p `x` p Hessian matrix at `vP`

returns
  1: successful, 0: function evaluation failed

Return value See `MaxBFGS`.
See also
  `MaxBFGS`, `GetMaxControl`, `MaxConvergenceMsg`, `Num1Derivative`, `Num2Derivative`

### MaxSimplex

    #import <maximize>
    MaxSimplex(const func, const avP, const adFunc, vDelta);

`func`
  in: a function computing the function value
`avP`
  in: address of p `x` 1 matrix with starting values
  out: p `x` 1 matrix with coefficients at convergence
`adFunc`
  in: address
  out: double, function value at convergence
`vDelta`
  in: 0, or the initial simplex

The supplied `func` argument should have the same format as in `MaxBFGS`.

Return value
  Returns the status of the iterative process, as documented under `MaxBFGS`.

### Num1Derivative, Num2Derivative

    #import <maximize>
    Num1Derivative(const func, vP, const avScore);
    Num2Derivative(const func, vP, const amHessian);

`func`
  in: a function computing the function value, optionally with derivatives
`vP`
  in: p `x` 1 matrix with parameter values
`mHessian`
  in: p `x` p matrix, initial Hessian
`avScore`
  in: an address
  out: p `x` 1 matrix with 1st derivatives at `vP`
`amHessian`
  in: an address
  out: p `x` p matrix with 2nd derivatives at `vP`

The supplied `func` argument should have the format as documented under `MaxBFGS`.

Return value
  Returns 1 if successful, 0 otherwise.
See also
  `MaxBFGS`

### NumJacobian

    #import <maximize>
    NumJacobian(const func, vU, const amJacobian);

`func`
  in: function mapping from restricted to unrestricted parameters
`vU`
  in: of u `x` 1 matrix with parameters
`amJacobian`
  in: address
  out: r `x` u Jacobian matrix corresponding to mapping

The supplied `func` argument should have the following format:

    func(const avR, const vU);

`avR`
  in: address
  out: r `x` 1 matrix with restricted coefficients
`vU`
  in: u `x` 1 matrix with unrestricted coefficients returns 1: successful, 0: function evaluation failed

Return value
  Returns 1 if successful, 0 otherwise.
See also
  `Num1Derivative`

### SolveQP

    #include <solveqp.h>
    SolveQP(const mG, const vG, const mA, const vB, const mC,
        const vD, const vLo, const vHi);
    SolveQPE(const mG, const vG, const mC, const vD);

`mG`
  in: n `x` n matrix G with quadratic weights, or n `x` 1 vector with diagonal of G
`vG`
  in: n `x` 1 vector g with linear weights
`mA`
  in: m `x` n matrix A with linear inequality constraints Ax>=b (may be empty)
`vB`
  in: m `x` 1 vector b with right-hand side for linear inequality constraints (empty if A is empty)
`mC`
  in: me `x` n matrix C with linear equality constraints cx=d (may be empty)
`vD`
  in: me `x` 1 vector d with right-hand side for linear equality constraints (empty if C is empty)
`vLo`
  in: n `x` 1 vector with lower bounds (may be empty)
`vHi`
  in: n `x` 1 vector with upper bounds (may be empty)

Return value
  `SolveQP` returns an array with three elements:
  - `[0]` integer return value:
    0: success
    1: initial point not feasible (should only be possible when `SolveQPIF` is called directly)
    2: maximum number of iterations reached
  - `[1]` n `x` 1 vector with solution x
  - `[2]` m* `x` 1 vector with Lagrange multipliers
  `SolveQPE` returns an array with three elements:
  - `[0]` n `x` 1 vector with solution x
  - `[1]` me `x` 1 vector with Lagrange multipliers
  - `[2]` p `x` 1 vector with index of redundant constraint (p=0 if all constraints were used)
Description
  `SolveQP` solves the quadratic program
  min x'Gx/2 + x'g, subject to:
  Ax >= b,
  Cx = d,
  xlo <= x <= xhi.
  using an active set method. The equality-constrained sub-problem is solved by `SolveQPE` using a null-space method.

  See `samples/maximize/solveqp1.ox` for an example.

## Probability package

### densbeta, densgamma, densgh, densgig, denskernel, densmises, denspoisson,

    #include <oxprob.h>
    densbeta(const ma, const a, const b);
    densgh(const x, const nu, const delta, const gamma, const beta);
    densgig(const x, const nu, const delta, const gamma);
    densgamma(const ma, const dr, const da);
    denskernel(const ma, const itype);
    densmises(const ma, const mu, const kappa);
    denspoisson(const ma, const mu);

`ma`
  in: arithmetic type
`a`,`b`
  in: arithmetic type, arguments for Beta distribution
`dr`
  in: arithmetic type
`da`
  in: arithmetic type
`mu`
  in: arithmetic type, von Mises: mean direction; Poisson: mean
`kappa`
  in: arithmetic type, dispersion parameter

Return value
  Returns the requested density at `ma` (the returned densities are between positive):
  - `densbeta`: Beta(a,b) density,
  - `densgamma`: Gamma(r, a) density,
  - `densgh`: GH(nu,delta,gamma,beta) density,
  - `densgig`: GIG(nu,delta,gamma) density,
  - `denskernel`: kernel densities, based on itype argument:
    - 'e' Epanechnikov (non-0 for x between -1 and 1)
    - 'b' Biweight (Quartic) (non-0 for x between -1 and 1)
    - 't' Triangular (non-0 for x between -1 and 1)
    - 'g' Gaussian (Normal)
    - 'r' Rectangular (Uniform) (non-0 for x between -1 and 1)
  - `densmises`: VM(mu, kappa) density (von Mises),
  - `denspoisson`: Poisson density,
  The return type is derived as follows:
  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest is scalar;
  - m `x` n matrix, when `ma` is a scalar, and the rest are m `x` n matrices;
  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest are m `x` n matrices;
  - double, when `ma` is scalar, and the rest is also scalar (int for `denst`).
See also
  `prob...`, `quan...`, `tail...`

### probbeta, probbvn, probgamma, probmises, probmvn, probpoisson,

    #include <oxprob.h>
    probbeta(const ma, const a, const b);
    probbvn(const da, const db, const drho);
    probgamma(const ma, const dr, const da);
    probmises(const ma, const mu, const kappa);
    probmvn(const mx, const msigma);
    probpoisson(const ma, const mu);

`ma`
  in: arithmetic type
`a`,`b`
  in: arithmetic type, arguments for Beta distribution
`dr`, `da`
  in: arithmetic type
`idf`
  in: int, degrees of freedom
`mu`
  in: arithmetic type, von Mises: mean location (use `M_PI` for symmetric between 0 and 2pi); Poisson: mean
`kappa`
  in: arithmetic type, dispersion parameter
`da`, `db`
  in: arithmetic type, upper limits of integration
`drho`
  in: arithmetic type, correlation coefficient
`mx`
  in: m `x` n matrix for n-variate normal, with m upper limits
`msigma`
  in: n `x` n variance matrix

Return value
  Returns the requested probabilities at `ma` (the returned probabilities are between zero and one):
  - `probbeta`: probabilities from Beta(a,b) distribution,
  - `probbvn`: probabilities from the bivariate normal distribution,
  - `probgamma`: probabilities from the Gamma distribution,
  - `probmises`: probabilities from the VM(mu, kappa) distribution (von Mises),
  - `probmvn`: probabilities from the multivariate normal distribution (currently up to trivariate),
  - `probpoisson`: cumulative probabilities from the Poisson distribution,

  The probabilities are accurate to about 10 digits. The return type for `probbeta`, `probgamma`, `probmises`, `probpoisson` is derived as follows:

  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest is scalar;
  - m `x` n matrix, when `ma` is a scalar, and the rest are m `x` n matrices;
  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest are m `x` n matrices;
  - double, when `ma` is scalar, and the rest is also scalar (int for `probt`).

  The return type for `probbvn` is a double if all arguments are scalar, or an m `x` n matrix if one or more arguments are an m `x` n matrix.

  The return type for `probmvn` is a double if m=1, or an m `x` 1 vector if m>1, where m is the number of rows of the first argument.

  Note that the von Mises distribution runs from 0 to 2pi, with mean direction mu (i.e. what tends to be written as VM(0,k) is VM(pi,k) here).

See also
  `bessel`, `betafunc`, `gammafunc`, `dens...`, `quan...`, `tail...`

### quanbeta, quangamma, quanmises,

    #include <oxprob.h>
    quanbeta(const ma, const a, const b);
    quangamma(const ma, const dr, const da);
    quanmises(const ma, const mu, const kappa);

`ma`
  in: arithmetic type, probabilities: all values must be between 0 and 1
`a`,`b`
  in: arithmetic type, arguments for Beta distribution
`dr`, `da`
  in: arithmetic type
`idf`
  in: int, degrees of freedom
`kappa`
  in: arithmetic type, dispersion parameter

Return value
  Returns the requested quantiles (inverse probability function; percentage points) at `ma`:
  - `quanbeta`: quantiles from Beta(a,b) distribution,
  - `quangamma`: quantiles from the Gamma distribution,
  - `quanmises`: quantiles from the VM(mu, kappa) distribution (von Mises),

  The quantiles are accurate to about 10 digits. The return type is derived as follows:

  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest is scalar;
  - m `x` n matrix, when `ma` is a scalar, and the rest are m `x` n matrices;
  - m `x` n matrix, when `ma` is an m `x` n matrix, and the rest are m `x` n matrices;
  - double, when `ma` is scalar, and the rest is also scalar.
See also
  `dens...`, `prob...`, `tail...`

### ranbeta, ranbinomial, ranbrownianmotion, ranchi, randirichlet, ranexp, ranf, rangamma, rangh, rangig, raninvgaussian, ranlogarithmic, ranlogistic, ranlogn, ranmises, ranmultinomial, rannegbin, ranpoisson, ranpoissonprocess, ranstable, ransubsample, ranshuffle, rant, ranuorder, ranwishart

    #include <oxprob.h>
    ranbeta(const r, const c, const a, const b);
    ranbinomial(const r, const c, const n, const p);
    ranbrownianmotion(const r, const times);
    ranchi(const r, const c, const df);
    randirichlet(const r, const valpha);
    ranexp(const r, const c, const lambda);
    ranf(const r, const c, const df1, const df2);
    rangamma(const r, const c, const dr, const da);
    rangh(const r, const c, const nu, const delta, const gamma, const beta);
    rangig(const r, const c, const nu, const delta, const gamma);
    raninvgaussian(const r, const c, const mu, const lambda);
    ranlogarithmic(const r, const c, const alpha);
    ranlogistic(const r, const c);
    ranlogn(const r, const c);
    ranmises(const r, const c, const kappa);
    ranmultinomial(const n, const vp);
    rannegbin(const r, const c, const n, const p);
    ranpoisson(const r, const c, const mu);
    ranpoissonprocess(const r, const times, const mu);
    ranshuffle(const c, const x);
    ranstable(const r, const c, const alpha, const beta);
    ransubsample(const c, const n);
    rant(const r, const c, const df);
    ranuorder(const c);
    ranwishart(const r);

`r`
  in: int, number of rows
`c`
  in: int, number of columns
`a`,`b`
  in: arithmetic type, arguments for Beta distribution
`n`
  in: int, number of trials
`p`
  in: double, probability of success
`vp`
  in: column or row vector with c probabilities of success (must sum to one)
`lambda`
  in: double
`df`
  in: double, degrees of freedom
`df1`
  in: double, degrees of freedom in the numerator
`df2`
  in: double, degrees of freedom in the denominator
`dr`
  in: double
`da`
  in: double
`mu`
  in: double, mean
`kappa`
  in: arithmetic type, dispersion parameter (mean direction is pi)
`alpha`
  in: double
`beta`
  in: double
`valpha`
  in: vector with c+1 shape parameters for Dirichlet distribution
`times`
  in: column or row vector with c time points (must be non-decreasing)
`x`
  in: column or row vector to sample from

Return value
  The following return a `r` by `c` matrix of random numbers from the selected distribution:
  - `ranbeta` : Beta(a,b) distribution,
  - `ranbinomial`: Binomial(n,p) distribution,
  - `ranbrownianmotion` : r realizations of the Brownian motion,
  - rancauchy: equals `rant(r, c, 1)`,
  - `ranchi` : chi2(df) distribution,
  - `randirichlet` : Dirichlet(\alpha1,...,\alphac+1) distributed random numbers (each row is a realization of the c-variate random variable),
  - `ranexp` : exp(\lambda) distribution with mean 1 / \lambda,
  - `ranf` : F(df1, df2) distribution,
  - `rangamma` : Gamma(r,a) distribution,
  - `rangh`: GH(nu,delta,gamma,beta) distribution,
  - `rangig`: GIG(nu,delta,gamma) distribution,
  - `raninvgaussian` : Inverse Gaussian(mu,lambda) distribution,
  - `ranlogarithmic` : logarithmic distribution,
  - `ranlogistic` : logistic distribution,
  - `ranlogn` : log normal distribution,
  - `ranmises`: VM(pi, kappa) distribution (von Mises),
  - `rannegbin`: Negative binomial(n,p) distribution,
  - `ranpoisson` : Poisson(mu) distribution,
  - `ranpoissonprocess` : r realizations of the Poisson process,
  - `ranstable` : S(alpha,beta) distribution,
  - `rant` : Student t(df) distribution, degrees of freedom need not be integer,
  - `ranwishart` : Wishart(1, Ir) distributed random drawing, returns an r `x` r matrix.

  The matrix is filled by row. Note that, if both `r` and `c` are 1, the return value is a scalar of type double!

  The following return a `1` by `c` matrix of random numbers from the selected distribution:

  - `ranmultinomial`: Multinomial(n,p1,...,pc) distribution,
  - `ranshuffle`: draws c elements from x without replacement,
  - `ransubsample`: draws c numbers from the integers 0,...,n-1 without replacement,
  - `ranuorder`: generates c uniform order statistics.

  The von Mises is generated between 0 and 2pi, with mean direction pi, corresponding to `probmises(.,M_PI,kappa)`. To use a different mean:

          y = fmod(ranmises(r, c, kappa) + mu, M_2PI);

  `M_PI` and `M_2PI` require `oxfloat.h`.
See also
  `rann`, `ranseed`, `ranu`

## QuadPack

QuadPack is a Fortran library for univariate numerical integration (`quadrature') using adaptive rules. Five functions are exported to Ox from a dynamic link library, using the header file `quadpack.h`.

### QAG, QNG, QAGS, QAGP, QAGI, QPEPS,

    #include <quadpack.h>
    QNG (const func, const a, const b, const aresult,
        const aabserr);
    QAG (const func, const a, const b, const key,
        const aresult, const aabserr);
    QAGS(const func, const a, const b, const aresult,
        const aabserr);
    QAGP(const func, const a, const b, const vpoints,
        const aresult, const aabserr);
    QAGI(const func, const bound, const inf, const aresult,
        const aabserr);
    QPEPS(const epsabs, const epsrel);

`func`
  in: function to integrate; func must be a function of one argument (a double), returning a double
`a`
  in: double, lower limit of integration
`b`
  in: double, upper limit of integration
`key`
  in: int, key for choice of local integration rule, which determines the number of points in the Gauss-Kronrod pair: 1 or less (7-15 points), 2 (10-21 points), 3 (15-31 points), 4 (20-41 points), 5 (25-51 points), 6 or more (30-61 points).
`vpoints`
  in: row vector with singularities of integrand
`bound`
  in: double, lower bound (`key` == 1) or upper bound (`key` == -1)
`inf`
  in: int, 1: integral from b to infinity, -1: integral from minus infinity to b, 2: integral from minus infinity to infinity
`aresult`
  in: address of variable
  out: double, approximation to the integral
`aabserr`
  in: address of variable
  out: double, estimate of the modulus of the absolute error
`epsabs`
  in: double, absolute accuracy requested (the default value is 0)
`epsrel`
  in: double, relative accuracy requested (the default value is 1e-10)

Return value
  Result of the QuadPack routine:
  - 0: normal and reliable termination of routine;
  - 1: maximum number of steps has been executed;
  - 2: roundoff error prevents reaching the desired tolerance;
  - 3: extremely bad integrand behaviour prevents reaching tolerance;
  - 4: algorithm does not converge;
  - 5: integral is probably convergent or slowly divergent;
  - 6: invalid input;
  - 10: not enough memory;
Description
  `QNG`: simple non-adaptive automatic integrator for a smooth integrand.
  `QAG`: simple globally adaptive Gauss-Kronrod-based integrator, with choice of formulae.
  `QAGS`: globally adaptive integrator with extrapolation, which can handle integrand singularities of several types.
  `QAGP`: as `QAGS`, but allows the user to specify singularities, discontinuities and other difficulties of the integrand.
  `QAGI`: as `QAGS`, but handles integration over infinite integrals.

  Full documentation is in Piessens et al. (1983), QUADPACK, A Subroutine Package for Automatic Integration, Springer-Verlag.

Ox version 3.00. © JA Doornik This file last changed
.