# Ox Standard Packages Reference

## Chapter contents:

| Arma |
arma0 | armaforc | armagen | armavar
diffpow | modelforc | pacf
| Maximization |
GetMaxControl | GetMaxControlEps | MaxBFGS | MaxControl
MaxControlEps | MaxConvergenceMsg | MaxSQPF | MaxNewton
MaxSimplex | MaxSQP | MaxSQPF | Num1Derivative
Num2Derivative | NumJacobian | SolveNLE | SolveQP
| Probability |
densbeta | densbinomial | denscauchy | densexp
densextremevalue | densgamma | densgeometric | densgh
densgig | denshypergeometric | densinvgaussian | denskernel
denslogarithmic | denslogistic | denslogn | densmises
densnegbin | denspareto | denspoisson | densweibull
probbeta | probbinomial | probbvn | probcauchy
probexp | probextremevalue | probgamma | probgeometric
probhypergeometric | probinvgaussian | problogarithmic | problogistic
problogn | probmises | probmvn | probnegbin
probpareto | probpoisson | probweibull
quanbeta | quanbinomial | quancauchy | quanexp
quanextremevalue | quangamma | quangeometric | quanhypergeometric
quaninvgaussian | quanlogarithmic | quanlogistic | quanlogn
quanmises | quannegbin | quanpareto | quanpoisson
quanweibull
ranbeta | ranbinomial | ranbrownianmotion | rancauchy
ranchi | randirichlet | ranexp | ranextremevalue
ranf | rangamma | rangeometric | rangh
rangig | ranhypergeometric | ranlogarithmic | ranlogistic
ranlogn | ranmises | ranmultinomial | rannegbin
ranpareto | ranpoisson | ranpoissonprocess | ranshuffle
ranstable | ransubsample | rant
ranuorder | ranwishart | ranweibull
| QuadPack |
QNG | QAG | QAGI | QAGP
QAGS | QAWC | QAWF | QAWO
QAWS | QPEPS | QPWARN

## 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`

Examples
  See `ox/samples/maximize/`.

### 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`, `MaxNewton`, `MaxSQP`, `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`

### 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`

Examples
  See `ox/samples/maximize/`.

### 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 a p `x` 1 matrix with the initial simplex (if 0 is specified, the score is used for 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`.

### MaxSQP, MaxSQPF

    #import <maxsqp>
    MaxSQP(const func, const avP, const adFunc, const amHessian,
        const fNumDer, const cfunc_gt0, const cfunc_eq0, vLo, vHi);
    MaxSQP(const func, const avP, const adFunc, const amHessian,
        const fNumDer, const cfunc_gt0, const cfunc_eq0, vLo, vHi,
    	const cfunc_gt0_jac, const cfunc_eq0_jac, ...);
    MaxSQPF(const func, const avP, const adFunc, const amHessian,
        const fNumDer, const cfunc_gt0, const cfunc_eq0, vLo, vHi);
    MaxSQPF(const func, const avP, const adFunc, const amHessian,
        const fNumDer, const cfunc_gt0, const cfunc_eq0, vLo, vHi,
    	const cfunc_gt0_jac, const cfunc_eq0_jac, ...);

`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 approximation matrix B
`fNumDer`
  in: 0: `func` provides analytical second derivatives 1: use numerical second derivatives
`vLo`
  in: p `x` 1 matrix with lower bounds, or <> if there are no lower bounds
`vHi`
  in: p `x` 1 matrix with upper bounds, or <> if there are no lower bounds

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

The `cfunc_gt0` argument can be zero, or a function evaluating the nonlinear constraints (which will be constrained to be positive) with the following format:

    cfunc_gt0(const avF, const vP);

`avF`
  in: address
  out: m `x` 1 matrix with constraints at `vP`
`vP`
  in: p `x` 1 matrix with coefficients

returns
  1: successful, 0: constraint evaluation failed

The `cfunc_eq0` argument can be zero, or a function evaluating the nonlinear constraints (which will be constrained to zero) with the following format:

    cfunc_eq0(const avF, const vP);

`avF`
  in: address
  out: me `x` 1 matrix with equality constraints at `vP`
`vP`
  in: p `x` 1 matrix with coefficients

returns
  1: successful, 0: constraint evaluation failed

The `cfunc_gt0_jac` and `cfunc_eq0_jac` are optional functions that return the analytical Jacobian matrix of the constraints. They have the same format, returning in `avF` an m `x` p and an me `x` p matrix respectively.

Return value See `MaxBFGS`.

Description `MaxSQP` implements a sequential quadratic programming technique to maximize a non-linear function subject to non-linear constraints, similar to Algorithm 18.7 in Nocedal and Wright (1999, Numerical Optimization).

`MaxSQPF` enforces all iterates to be feasible, using the Algorithm by Lawrence and Tits (2001, SIAM J. Optim.). The current version does not support equality constraints. If a starting point is infeasible, `MaxSQPF` will try to minimize the squared constraint violations to find a feasible point.

See also
  `MaxBFGS`, `GetMaxControl`, `MaxConvergenceMsg`, `Num1Derivative`, `Num2Derivative`

Examples
  See `ox/samples/maximize/`.

### 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`
Examples
  See `ox/samples/maximize/numder.ox`.

### 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`
Examples
  See `ox/samples/maximize/jacobian.ox`.

### SolveNLE

    #import <solvenle>
    SolveNLE(const func, const avX);
    SolveNLE(const func, const avX, iMode, funcJac, dEps1, dEps2,
        mxIter, iPrint, mxItInner);

`func`
  in: Ox function evaluating the nonlinear equations (see below)
`avX`
  in: address of n `x` 1 matrix with starting values
  out: n `x` 1 matrix with final coefficients
`iMode`
  in: int, mode of operation:

  -1 (default) | mode 1 if n < 80, else mode 3
  0 | tensor-Newton method using analytical Jacobian
  1 | tensor-Newton method using numerical Jacobian
  2 | tensor method using Broyden's approximation to Jacobian
  3 | large scale problem (tensor-gmres-Newton method, avoiding n `x` n Jacobian matrix)

`dEps1`
  in: double, eps1, default is 1e-4, use ≤ 0 to leave the current value unchanged (can also be set with `MaxControlEps`)
`dEps2`
  in: double, eps2, default is 5e-3, use ≤ 0 to leave the current value unchanged (can also be set with `MaxControlEps`)
`mxIter`
  in: int, maximum number of iterations; default is 1000, use -1 to leave the current value unchanged (can also be set with `MaxControls`)
`iPrint`
  in: int, print results every iPrint'th iteration; default is 0, use -1 to leave the current value unchanged (can also be set with `MaxControls`)
`mxItInner`
  in: int, number of inner iterations for large scale problems, default is max(50, 10 * log10(n))

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

    func(const avF, const vX)

`avF`
  in: address
  out: n `x` 1 matrix with nonlinear system f(x) evaluated at x
`vX`
  in: n `x` 1 matrix with x values

returns
  1: successful, 0: function evaluation failed

When the analytical Jacobian is used, the `funcJac` argument should have the following format:

    funcJac(const amJac, const vX)

`amJac`
  in: address
  out: n `x` n matrix with Jacobian matrix evaluated at `vX`
`vX`
  in: n `x` 1 matrix with X values unrestricted coefficients returns 1: successful, 0: function evaluation failed

returns
  1: successful, 0: function evaluation failed

Return value
  Returns the status of the iterative process:
  `MAX_CONV`
    Strong convergence norm(f(x)) < 0.001eps1.
  `MAX_WEAK_CONV`
    Weak convergence (no improvement in line search) The step length has become too small and norm(f(x)) < eps2.
  `MAX_MAXIT`
    No convergence (maximum no of iterations reached)
  `MAX_LINE_FAIL`
    No convergence (no improvement in line search) The step length has become too small and weak convergence was not achieved.
  `MAX_FUNC_FAIL`
    No convergence (function evaluation failed)
  `MAX_NOCONV`
    No convergence Probably not yet attempted to solve the system.

Description
  Solves a system f(x) of n nonlinear equations in n unknowns.

See also
  `MaxControls`, `MaxControlsEps`

Examples
  See `ox/samples/maximize/solvenle1.ox` and `ox/samples/maximize/solvenle2.ox`.

### SolveQP

    #import <solveqp>
    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);
    SolveQPS(const sFile, const iVerbose)
    SolveQPS(const sFile, const iVerbose, const fnSolveQP)

`mG`
  in: n `x` n symmetric 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)
`sFile`
  in: string with `.qps` file name
`iVerbose`
  in: int, 0 for no output, 1 for one line summary output, 2 to print all matrices and results.
`fnSolveQP`
  in: (optional argument) QP solver with call syntax as `SolveQP`. If absent `SolveQP` is used.

Return value
  `SolveQP` returns an array with three elements:
  - `[0]` integer return value:
    0: success
    1: infeasible constraints
    2: maximum number of iterations reached
    3: Choleski decomposition failed
  - `[1]` n `x` 1 vector with solution x
  - `[2]` m* `x` 1 vector with Lagrange multipliers, m*=me+m+2n
    in order: equality constraints, inequality constraints, lower bounds, upper bounds.
  `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)
  `SolveQPS` returns an array with four elements: the first three as `SolveQP`, the fourth is the value of the objective function f(x).
Description
  `SolveQP` solves the quadratic program
  min f(x)=x'Gx/2 + x'g, subject to:
  Ax >= b,
  Cx = d,
  xlo <= x <= xhi.
  using an active set method. If G is not positive definite, a small number is added to its diagonal.

  `SolveQPE` solves an equality-constrained problem using a null-space method.

  `SolveQPS` solves a QP problem that is stored in a `.qps` file.

Examples
  See `samples/maximize/solveqp1.ox`.

## Probability package

### densbeta, densbinomial, denscauchy, densexp, densextremevalue, densgamma, densgeometric, densgh, densgig, denshypergeometric, densinvgaussian, denskernel, denslogarithmic, denslogistic, denslogn, densmises, densnegbin, denspareto, denspoisson, densweibull,

    #include <oxprob.h>
    densbeta(const ma, const a, const b);
    densbinomial(const ma, const n, const p);
    denscauchy(const ma);
    densexp(const ma, const lambda);
    densextremevalue(const ma, const alpha, const beta);
    densgamma(const ma, const dr, const da);
    densgeometric(const ma, const p);
    densgh(const ma, const nu, const delta, const gamma, const beta);
    densgig(const ma, const nu, const delta, const gamma);
    denshypergeometric(const ma, const n, const k, const m);
    densinvgaussian(const ma, const mu, const lambda);
    denskernel(const ma, const itype);
    denslogarithmic(const ma, const alpha);
    denslogistic(const ma, const alpha, const beta);
    denslogn(const ma);
    densmises(const ma, const mu, const kappa);
    densnegbin(const ma, const k, const p);
    denspareto(const ma, const k, const a);
    denspoisson(const ma, const mu);
    densweibull(const ma, const a, const b);

`ma`
  in: arithmetic type
`a`,`b`
  in: arithmetic type, arguments for Beta distribution
`dr`
  in: arithmetic type
`da`
  in: arithmetic type
`alpha`,`beta`
  in: arithmetic type, location and scale parameter
`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):

  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, probbinomial, probbvn, probcauchy, probexp, probextremevalue, probgamma, probgeometric, probhypergeometric, probinvgaussian, problogarithmic, problogistic, problogn, probmises, probmvn, probnegbin, probpareto, probpoisson, probweibull,

    #include <oxprob.h>
    probbeta(const ma, const a, const b);
    probbinomial(const ma, const n, const p);
    probbvn(const da, const db, const drho);
    probcauchy(const ma);
    probexp(const ma, const lambda);
    probextremevalue(const ma, const alpha, const beta);
    probgamma(const ma, const dr, const da);
    probgeometric(const ma, const p);
    probhypergeometric(const ma, const n, const k, const m);
    probinvgaussian(const ma, const mu, const lambda);
    problogarithmic(const ma, const alpha);
    problogistic(const ma, const alpha, const beta);
    problogn(const ma);
    probmises(const ma, const mu, const kappa);
    probmvn(const mx, const msigma);
    probnegbin(const ma, const k, const p);
    probpareto(const ma, const k, const a);
    probpoisson(const ma, const mu);
    probweibull(const ma, const a, const b);

`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
`alpha`,`beta`
  in: arithmetic type, location and scale 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 cumulative distribution functions at `ma` (P[X≤= x]; the returned probabilities are between zero and one):
  - `probbeta`: probabilities from Beta(a,b) distribution,
  - `probbinomial`: probabilities from Bin(n,p) distribution,
  - `probbvn`: probabilities from the bivariate normal distribution,
  - `probcauchy`: probabilities from the Cauchy distribution,
  - `probexp`: probabilities from the exp(lambda) distribution with mean 1 / lambda,
  - `probextremevalue`: probabilities from the Extreme Value (type I or Gumbel) distribution,
  - `probgamma`: probabilities from the Gamma distribution,
  - `probgeometric`: probabilities from the Geometric distribution,
  - `probhypergeometric`: probabilities from the Hypergeometric distribution,
  - `probinvgaussian`: probabilities from the Inverse Gaussian distribution,
  - `problogarithmic`: probabilities from the Logarithmic distribution,
  - `problogistic`: probabilities from the Logistic distribution,
  - `problogn`: probabilities from the Lognormal distribution,
  - `probmises`: probabilities from the VM(mu, kappa) distribution (von Mises),
  - `probmvn`: probabilities from the multivariate normal distribution (currently up to trivariate),
  - `probnegbin`: probabilities from the Negative Binomial distribution,
  - `probpareto`: probabilities from the Pareto distribution,
  - `probpoisson`: probabilities from the Poisson mu distribution,
  - `probweibull`: probabilities from the Weibull distribution.

  The functional forms are listed under the density functions.

  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, quanbinomial, quancauchy, quanexp, quanextremevalue, quangamma, quangeometric, quanhypergeometric, quaninvgaussian, quanlogarithmic, quanlogistic, quanlogn, quanmises, quannegbin, quanpareto, quanpoisson, quanweibull,

    #include <oxprob.h>
    quanbeta(const ma, const a, const b);
    quanbinomial(const ma, const n, const p);
    quancauchy(const ma);
    quanexp(const ma, const lambda);
    quanextremevalue(const ma, const alpha, const beta);
    quangamma(const ma, const dr, const da);
    quangeometric(const ma, const p);
    quanhypergeometric(const ma, const n, const k, const m);
    quaninvgaussian(const ma, const mu, const lambda);
    quanlogarithmic(const ma, const alpha);
    quanlogistic(const ma, const alpha, const beta);
    quanlogn(const mx);
    quanmises(const mp, const mu, const kappa);
    quannegbin(const ma, const k, const p);
    quanpareto(const ma, const k, const a);
    quanpoisson(const ma, const mu);
    quanweibull(const ma, const a, const b);

`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
`alpha`,`beta`
  in: arithmetic type, location and scale parameter
`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,
  - `quanbinomial`: quantiles from Bin(n,p) distribution,
  - `quancauchy`: quantiles from the Cauchy distribution,
  - `quanexp`: quantiles from the exp(lambda) distribution with mean 1 / lambda,
  - `quanextremevalue`: quantiles from the Extreme Value (type I or Gumbel) distribution,
  - `quangamma`: quantiles from the Gamma distribution,
  - `quangeometric`: quantiles from the Geometric distribution,
  - `quanhypergeometric`: quantiles from the Hypergeometric distribution,
  - `quaninvgaussian`: quantiles from the Inverse Gaussian distribution,
  - `quanlogarithmic`: quantiles from the Logarithmic distribution,
  - `quanlogistic`: quantiles from the Logistic distribution,
  - `quanlogn`: quantiles from the Lognormal distribution,
  - `quanmises`: quantiles from the VM(mu, kappa) distribution (von Mises),
  - `quannegbin`: quantiles from the Negative Binomial distribution,
  - `quanpareto`: quantiles from the Pareto distribution,
  - `quanpoisson`: quantiles from the Poisson mu distribution,
  - `quanweibull`: quantiles from the Weibull distribution.

  The functional forms are listed under the density functions.

  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, rancauchy, ranchi, randirichlet, ranexp, ranextremevalue, ranf, rangamma, rangeometric, rangh, rangig, ranhypergeometric, raninvgaussian, ranlogarithmic, ranlogistic, ranlogn, ranmises, ranmultinomial, rannegbin, ranpareto, ranpoisson, ranpoissonprocess, ranshuffle, ranstable, ransubsample, rant, ranuorder, ranwishart, ranweibull

    #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);
    rancauchy(const r, const c);
    ranchi(const r, const c, const df);
    randirichlet(const r, const valpha);
    ranexp(const r, const c, const lambda);
    ranextremevalue(const r, const c, const alpha, const beta);
    ranf(const r, const c, const df1, const df2);
    rangamma(const r, const c, const dr, const da);
    rangeometric(const r, const c, const p);
    rangh(const r, const c, const nu, const delta, const gamma,
        const beta);
    ranhypergeometric(const r, const c, const n, const k, const m);
    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);
    ranpareto(const r, const c, const k, const a);
    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);
    ranweibull(const r, const c, const a, const b);
    ranwishart(const r);

`r`
  in: int, number of rows
`c`
  in: int, number of columns
`a`,`b`
  in: double or r `x` c matrix, arguments for Beta distribution
`n`
  in: int, number of trials
`p`
  in: double, probability of success (rangeometric also allows or r `x` c matrix)
`vp`
  in: column or row vector with c probabilities of success (must sum to one)
`lambda`
  in: double or r `x` c matrix
`df`
  in: double or r `x` c matrix, degrees of freedom
`df1`
  in: double or r `x` c matrix, degrees of freedom in the numerator
`df2`
  in: double or r `x` c matrix, degrees of freedom in the denominator
`dr`
  in: double or r `x` c matrix
`da`
  in: double or r `x` c matrix
`mu`
  in: double or r `x` c matrix, mean
`kappa`
  in: double or r `x` c matrix, dispersion parameter (mean direction is pi)
`alpha`
  in: double or r `x` c matrix
`beta`
  in: double or r `x` c matrix
`nu`
  in: double, parameter for GH and GIG distributions
`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, with time steps specified by the vector `times`,
  - 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,
  - `ranextremevalue` : Extreme Value (type I or Gumbel) distribution
  - `ranf` : F(df1, df2) distribution,
  - `rangamma` : Gamma(r,a) distribution,
  - `rangeometric`: Geometric distribution,
  - `rangh`: GH(nu,delta,gamma,beta) distribution,
  - `rangig`: GIG(nu,delta,gamma) distribution,
  - `ranhypergeometric`: Hypergeometric 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,
  - `ranpareto`: Pareto(k,a) distribution,
  - `ranpoisson` : Poisson(mu) distribution,
  - `ranpoissonprocess` : r realizations of the Poisson process, with time steps specified by the vector `times`,
  - `ranstable` : S(alpha,beta) distribution,
  - `rant` : Student t(df) distribution, degrees of freedom need not be integer,
  - `ranweibull`: Weibull distribution.
  - `ranwishart` : Wishart(1, Ir) distributed random drawing, returns an r `x` r matrix.

  The functional forms are listed under the density functions.

  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. The main driver functions are exported to Ox from a dynamic link library, using the header file `quadpack.h`.

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

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

    #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);

`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 (`inf` == 1) or upper bound (`inf` == -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

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;
  An error message greater than 0 is reported unless switched off with `QPWARN`.
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.

### QAWO, QAWF, QAWS, QAWC

    #include <quadpack.h>
    QAWO(const func, const a, const b, const omega, const fcos,
    	const maxp1, const aresult, const aabserr);
    QAWF(const func, const a, const omega, const fcos, const limlst,
    	const maxp1, const aresult, const aabserr);
    QAWS(const func, const a, const b, const alpha, const beta,
    	const type, const aresult, const aabserr);
    QAWC(const func, const a, const b, const c, const aresult,
    	const aabserr);

`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
`omega`
  in: double, factor in cosine or sine function
`fcos`
  in: int, 1: function to integrate is cos(omega*x)*f(x), else it is sin(omega*x)*f(x)
`maxp1`
  in: int, upper bound on the number of Chebyshev moments which can be stored.
`limlst`
  in: int, upper bound on the number of cycles (must be >= 3).
`alpha,beta`
  in: double, powers in w(x), both >-1.
`itype`
  in: int, 1: v(x)=1; 2: v(x)=log(x-a); 3: v(x)=log(b-x); 4: v(x)=log(x-a)*log(b-x).
`c`
  in: double, term for Cauchy principal value (!= a and != b).
`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

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;
  An error message greater than 0 is reported unless switched off with `QPWARN`.
Description
  `QAWO`: integrates cos(omega*x)*f(x) or sin(omega*x)*f(x) over a finite interval (a,b).
  `QAWF`: Fourier cosine or Fourier sine transform of f(x), from a to infinity. (`QAWF` returns error 6 if epsabs is zero, use `QPEPS` to change the value of epsabs.)
  `QAWS`: integrates w(x)*f(x) over a finite interval (a,b), where w(x) = [(x-a)^alpha]*[(b-x)^beta]*v(x), and v(x) depends on the `itype` argument.
  `QAWC`: Cauchy principal value of f(x)/(x-c) over a finite interval (a,b) and for user-determined c.

### QPEPS, QPWARN

    #include <quadpack.h>
    QPEPS(const epsabs, const epsrel);
    QPWARN(const ion);

`epsabs`
  in: double, absolute accuracy requested (the default value is 0)
`epsrel`
  in: double, relative accuracy requested (the default value is 1e-10)
`ion`
  in: int, 1: print warning and error messages (the default), or 0: don't print

No return
Description
  `QPEPS`: sets the accuracy for all QuadPack functions.
  `QPWARN`: controls whether warning/error messages are printed or not.

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