# Ox Standard Packages Reference

## Chapter contents:

| Arma |
arma0 | armaforc | armagen | armavar
diffpow | modelforc | pacf
| Maximization Control |
CMaxControl class | GetMaxControl | GetMaxControlEps | MaxControl
MaxControlEps | MaxConvergenceMsg
| Maximization Functions |
MaxBFGS | MaxNewton | MaxSimplex | MaxSQP
MaxSQPF | Num1Derivative | Num1Derivative_parallel | Num2Derivative
Num2Derivative_parallel | 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 | ranindex | 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

The Arma package implements functions which are commonly used in autoregressive-moving average models. The Arma package requires the header file `arma.h`. Note that the Arma package uses the convention of writing the AR and MA coefficients on the right-hand side with a positive sign.

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

### Maximization control

### CMaxControl Class

`CMaxControl(const iOptions=0);`
  Constructor. Only possible option is `CMaxControl::PARALLEL_SCORE`.
`GetControl();`
  returns `{ mxIter, iPrint, bCompact }`.
`GetEps();`
  returns `{ dEps1, dEps2 }`.
`GetIterationCount();`
  returns the iteration count.
`GetResult();`
  returns the convergence code.
`SetControl(const mxIter, const iPrint=-1, const bCompact=-1);`
  See `MaxControl`.
`SetEps(const dEps1, const dEps2=-1);`
  See `GetMaxControl`.
`SetIterationCount(const cIter);`
  Sets the iteration count.
`SetOptions(const iOptions);`
  Only possible option is `CMaxControl::PARALLEL_SCORE`.
`SetResult(const iResult);`
  Sets the convergence code.

The `CMaxControl` class manages the configurable maximization options in a more convenient way. The added flexibility is that parallel numerical scores can be used, and the number of iterations retrieved upon convergence.

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

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

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

### Maximization functions

### MaxBFGS

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

  `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
  `objMaxCtrl`
    in: CMaxControl object (optional argument)
    out: updated to reflect status and iteration count.
  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`, CMaxControl, `MaxConvergenceMsg`, `Num1Derivative`, `Num2Derivative`
Examples
  See `ox/samples/maximize/`.

### MaxNewton

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

  `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
  `objMaxCtrl`
    in: CMaxControl object (optional argument)
    out: updated to reflect status and iteration count.
  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`, CMaxControl, `MaxConvergenceMsg`, `Num1Derivative`, `Num2Derivative`
Examples
  See `ox/samples/maximize/`.

### MaxSimplex

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

  `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)
  `objMaxCtrl`
    in: CMaxControl object (optional argument)
    out: updated to reflect status and iteration count.

  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 first derivatives 1: use numerical first 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, Num1Derivative_parallel Num2Derivative, Num2Derivative_parallel

    #import <maximize>
    Num1Derivative(const func, vP, const avScore);
    Num1Derivative_parallel(const func, vP, const avScore);
    Num2Derivative(const func, vP, const amHessian);
    Num2Derivative_parallel(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, ranindex, 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);
    ranindex(const c);
    ranindex(const c, const n);
    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 n, const p);

  `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(n, Ip) distributed random drawing, returns an p `x` p 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:

  - `ranindex(c)`: draws c numbers from 0,...,c-1 without replacement,
  - `ranindex(c,n)`: draws c numbers from 0,...,n-1 without replacement (this is the same as `ranshuffle(c, range(0,n-1))`),
  - `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, returning the ordered numbers (the return value is sorted, so `ransubsample(n,n)` just returns 0,...,n-1),
  - `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 7.00. © JA Doornik This file last changed
.