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 \phi₁, \phi₂, ..., \phi_p followed by MA coefficients \theta₁, \theta₂, ..., \theta_q, 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:

e₀ = 0,

e₁ = a₁ - \phi₁ a₀ - \theta₁ e₀,

e₂ = a₂ - \phi₁ a₁ - \theta₁ e₁ - \theta ₁ e₀, 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 \phi₁, \phi₂, ..., \phi_p followed by MA coefficients \theta₁, \theta₂, ..., \theta_q, 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:

ahat₀ = x₀ + \phi₁ a_T-1 + \phi₂ a_T-2 + \theta₁ e_T-1 + \theta₂ e_T-2,

ahat₁ = x₁ + \phi₁ ahat₀ + \phi₂ a_T-1 + \theta₂ e_T-1,

ahat₂ = x₂ + \phi₁ ahat₁ + \phi₂ ahat₀, 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 \phi₁, \phi₂, ..., \phi_p followed by MA coefficients \theta₁, \theta₂, ..., \theta_q, 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:

a₀ = x₀,

a₁ = x₁ + \phi₁ a₀ + e₁ + \theta₁ e₀,

a₂ = x₂ + \phi₁ a₁ + e₂ + \theta₁ e₁ + \theta₂ e₀, 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 \phi₁, \phi₂, ..., \phi_p followed by MA coefficients \theta₁, \theta₂, ..., \theta_q, 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, \sigma² _\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-T₀) 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-T₀) 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 = eps₁.

MAX_WEAK_CONV

Weak convergence (no improvement in line search) The step length s_i has become too small. The convergence test ) was passed, using tolerance eps = eps₂.

MAX_MAXIT

No convergence (maximum no of iterations reached)

MAX_LINE_FAIL

No convergence (no improvement in line search) The step length s_i has become too small. The convergence test was not passed, using tolerance eps = eps₂.

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:

eps₁=1e-4, eps₂=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, eps₁, default is 1e-4, use ≤ 0 to leave the current value unchanged

dEps2

in: double, eps₂, 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

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

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 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: m_e 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 m_e 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 nxn Jacobian matrix)

dEps1

in: double, eps₁, default is 1e-4, use ≤ 0 to leave the current value unchanged (can also be set with MaxControlEps)

dEps2

in: double, eps₂, 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.001eps₁.

MAX_WEAK_CONV

Weak convergence (no improvement in line search) The step length has become too small and norm(f(x)) < eps₂.

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: m_e x n matrix C with linear equality constraints cx=d (may be empty)

vD

in: m_e 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*=m_e+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] m_e 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,
x_lo <= x <= x_hi. 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 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 : chi²(df) distribution,
• randirichlet : Dirichlet(\alpha₁,...,\alpha_c+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, I_p) 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:

• ranmultinomial: Multinomial(n,p₁,...,p_c) 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 (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 6.10. © JA Doornik This file last changed .