Ox Standard Packages Reference

Chapter contents:

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 or H x 1 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 or T x 1 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

FindZero

#import <maxscalar>
FindZero(const Func, const dX, dStep, const dTol, const mxIter);
Func
in: a function of one argument returns the function value
Func
in: function to find root of: Func(x) returns the function value
dX
in: double, center of initial bracket
dStep
in: double, determines bracket [dX + 2^i * dStep,dX + 2^i * dStep] with sign change (default 0.1)
dTol
in: convergence tolerance (default 1e-13)
mxIter
in: int, maximum no of iterations (default 1000)
Return value Returns the root using the Brent (1973) algorithm to find root of a function of one variable without using derivatives.
Examples
See ox/samples/maximize/test_maxscalar.ox.

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/.

MaxScalarBrent

MaxScalarPowell

#import <maxscalar>
MaxScalarBrent(const Func, const dP0, const dP1, dF0, dF1,
    const dA, const dB, const dTol, const mxIter)
MaxScalarPowell(const Func, const dP0, const dP1, dF0, dF1,
    const dA, const dB, const dTol, const mxIter)
Func
in: a function of one argument returns the function value
Func
in: function to find root of: Func(x) returns the function value
dP0
in: double, first starting value
dP1
in: double, second starting value
dF0
in: double, function value at dP0 (default .NaN to force evaluation)
dF1
in: double, function value at dP1 (default .NaN to force evaluation)
dA,dB
in: doubles, defining interval to maximize over (default -100,100). Must have dP0,dP1 in (dA,dB), but this is not checked!
dTol
in: convergence tolerance (default 1e-7)
mxIter
in: int, maximum no of iterations (default 1000)
Return value
Returns an array with four elements: {argument at maximum, function value at maximum, no of function calls, no of iterations}.
Examples
See ox/samples/maximize/test_maxscalar.ox.

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
0tensor-Newton method using analytical Jacobian
1tensor-Newton method using numerical Jacobian
2tensor method using Broyden's approximation to Jacobian
3large scale problem (tensor-gmres-Newton method, avoiding nxn 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):
./oxdens1.png
./oxdens2.png
./oxdens3.gif
./oxdens4.gif
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 rxc matrix, arguments for Beta distribution
n
in: int, number of trials
p
in: double, probability of success (rangeometric also allows or rxc matrix)
vp
in: column or row vector with c probabilities of success (must sum to one)
lambda
in: double or rxc matrix
df
in: double or rxc matrix, degrees of freedom
df1
in: double or rxc matrix, degrees of freedom in the numerator
df2
in: double or rxc matrix, degrees of freedom in the denominator
dr
in: double or rxc matrix
da
in: double or rxc matrix
mu
in: double or rxc matrix, mean
kappa
in: double or rxc matrix, dispersion parameter (mean direction is pi)
alpha
in: double or rxc matrix
beta
in: double or rxc 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 9.00. © JA Doornik This file last changed .