default.QP {monomvn} | R Documentation

## 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 S 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 t(w) * S * w,

subject to the constraints that all 0 <= w[i] <=1, for i=1:m,

sum(w) = 1,

and where S 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 - t(w) * mu + 0.5 * t(w) * S * w,

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

t(mu) * w >= 1,

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

-t(mu) * w >= -2,

i.e., with alternating sign on the linear part, so that each sample of t(mu)*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 S 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`
`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 bobby@statslab.cam.ac.uk

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

[Package monomvn version 1.8 Index]