default.QP {monomvn}

Generating a default Quadratic Program for bmonomvn

Description

This function generates a default “minimum variance” Quadratic Program in order to obtain samples of the solution under the posterior for parameters \mu and \Sigma obtained via bmonomvn. The list generated as output has entries similar to the inputs of solve.QP from the quadprog package

Usage

default.QP(m, dmu = FALSE, mu.constr = NULL)

Arguments

m

the dimension of the solution space; usually ncol(y) or equivalently length(mu), ncol(S) and nrow(S) in the usage of bmonomvn

dmu

a logical indicating whether dvec should be replaced with samples of \mu; see details below

mu.constr

a vector indicating linear constraints on the samples of \mu to be included in the default constraint set. See details below; the default of NULL indicates none

Details

When bmonomvn(y, QP=TRUE) is called, this function is used to generate a default Quadratic Program that samples from the argument w such that

\min_w w^\top \Sigma w,

subject to the constraints that all 0\leq w_i \leq 1, for i=1,\dots,m,

\sum_{i=1}^m w_i = 1,

and where \Sigma is sampled from its posterior distribution conditional on the data y. Alternatively, this function can be used as a skeleton to for adaptation to more general Quadratic Programs by adjusting the list that is returned, as described in the “value” section below.

Non-default settings of the arguments dmu and mu.constr augment the default Quadratic Program, described above, in two standard ways. Specifying dvec = TRUE causes the program objective to change to

\min_w - w^\top \mu + \frac{1}{2} w^\top \Sigma w,

with the same constraints as above. Setting mu.constr = 1, say, would augment the constraints to include

\mu^\top w \geq 1,

for samples of \mu from the posterior. Setting mu.constr = c(1,2) would augment the constraints still further with

-\mu^\top w \geq -2,

i.e., with alternating sign on the linear part, so that each sample of \mu^\top w must lie in the interval [1,2]. So whereas dmu = TRUE allows the mu samples to enter the objective in a standard way, mu.constr (!= NULL) allows it to enter the constraints.

The accompanying function monomvn.solve.QP can act as an interface between the constructed (default) QP object, and estimates of the covariance matrix \Sigma and mean vector \mu, that is identical to the one used on the posterior-sample version implemented in bmonomvn. The example below, and those in the documentation for bmonomvn, illustrate how this feature may be used to extract mean and MLE solutions to the constructed Quadratic Program

Value

This function returns a list that can be interpreted as specifying the following arguments to the solve.QP function in the quadprog package. See solve.QP for more information of the general specification of these arguments. In what follows we simply document the defaults provided by default.QP. Note that the Dmat argument is not, specified as bmonomvn will use samples from S (from the posterior) instead

m

length(dvec), etc.

dvec

a zero-vector rep(0, m), or a one-vector rep(1, m) when dmu = TRUE as the real dvec that will be used by solve.QP will then be dvec * mu

dmu

a copy of the dmu input argument

Amat

a matrix describing a linear transformation which, together with b0 and meq, describe the constraint that the components of the sampled solution(s), w, must be positive and sum to one

b0

a vector containing the (RHS of) in/equalities described by the these constraints

meq

an integer scalar indicating that the first meq constraints described by Amat and b0 are equality constraints; the rest are >=

mu.constr

a vector whose length is one greater than the input argument of the same name, providing bmonomvn with the number

mu.constr[1] = length(mu.constr[-1])

and location mu.constr[-1] of the columns of Amat which require multiplication by samples of mu

The $QP object that is returned from bmonomvn will have the following additional field

o

an integer vector of length m indicating the ordering of the rows of $Amat, and thus the rows of solutions $W that was used in the monotone factorization of the likelihood. This field appears only after bmonomvn returns a QP object checked by the internal function check.QP

Author(s)

Robert B. Gramacy [e-mail]

See Also

bmonomvn and solve.QP in the quadprog package, monomvn.solve.QP

Examples

## generate N=100 samples from a 10-d random MVN
## and randomly impose monotone missingness
xmuS <- randmvn(100, 20)
xmiss <- rmono(xmuS$x)

## set up the minimum-variance (default) Quadratic Program
## and sample from the posterior of the solution space
qp1 <- default.QP(ncol(xmiss))
obl1 <- bmonomvn(xmiss, QP=qp1)
bm1 <- monomvn.solve.QP(obl1$S, qp1) ## calculate mean
bm1er <- monomvn.solve.QP(obl1$S + obl1$mu.cov, qp1) ## use estimation risk
oml1 <- monomvn(xmiss)
mm1 <- monomvn.solve.QP(oml1$S, qp1) ## calculate MLE

## now obtain samples from the solution space of the
## mean-variance QP
qp2 <- default.QP(ncol(xmiss), dmu=TRUE)
obl2 <- bmonomvn(xmiss, QP=qp2)
bm2 <- monomvn.solve.QP(obl2$S, qp2, obl2$mu) ## calculate mean
bm2er <- monomvn.solve.QP(obl2$S + obl2$mu.cov, qp2, obl2$mu) ## use estimation risk
oml2 <- monomvn(xmiss)
mm2 <- monomvn.solve.QP(oml2$S, qp2, oml2$mu) ## calculate MLE

## now obtain samples from minimum variance solutions
## where the mean weighted (samples) are constrained to be
## greater one
qp3 <- default.QP(ncol(xmiss), mu.constr=1)
obl3 <- bmonomvn(xmiss, QP=qp3)
bm3 <- monomvn.solve.QP(obl3$S, qp3, obl3$mu) ## calculate mean
bm3er <- monomvn.solve.QP(obl3$S + obl3$mu.cov, qp3, obl3$mu) ## use estimation risk
oml3 <- monomvn(xmiss)
mm3 <- monomvn.solve.QP(oml3$S, qp3, oml2$mu) ## calculate MLE

## plot a comparison
par(mfrow=c(3,1))
plot(obl1, which="QP", xaxis="index", main="Minimum Variance")
points(bm1er, col=4, pch=17, cex=1.5) ## add estimation risk
points(bm1, col=3, pch=18, cex=1.5) ## add mean
points(mm1, col=5, pch=16, cex=1.5) ## add MLE
legend("topleft", c("MAP", "posterior mean", "ER", "MLE"), col=2:5,
       pch=c(21,18,17,16), cex=1.5)
plot(obl2, which="QP", xaxis="index", main="Mean Variance")
points(bm2er, col=4, pch=17, cex=1.5) ## add estimation risk
points(bm2, col=3, pch=18, cex=1.5) ## add mean
points(mm2, col=5, pch=16, cex=1.5) ## add MLE
plot(obl3, which="QP", xaxis="index", main="Minimum Variance, mean >= 1")
points(bm3er, col=4, pch=17, cex=1.5) ## add estimation risk
points(bm3, col=3, pch=18, cex=1.5) ## add mean
points(mm3, col=5, pch=16, cex=1.5) ## add MLE

## for a further comparison of samples of the QP solution
## w under Bayesian and non-Bayesian monomvn, see the
## examples in the bmonomvn help file