Description Usage Arguments Details Value Author(s) See Also Examples
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
1  default.QP(m, dmu = FALSE, mu.constr = NULL)

m 
the dimension of the solution space; usually

dmu 
a logical indicating whether 
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 
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.
Nondefault 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 posteriorsample 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
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 

dvec 
a zerovector 
dmu 
a copy of the 
Amat 
a 
b0 
a vector containing the (RHS of) in/equalities described by the these constraints 
meq 
an integer scalar indicating that the first 
mu.constr 
a vector whose length is one greater than the
input argument of the same name, providing
and location 
The $QP
object that is returned from bmonomvn
will have the following additional field
o 
an integer vector of length 
Robert B. Gramacy rbg@vt.edu
bmonomvn
and solve.QP
in the quadprog package, monomvn.solve.QP
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53  ## generate N=100 samples from a 10d random MVN
## and randomly impose monotone missingness
xmuS < randmvn(100, 20)
xmiss < rmono(xmuS$x)
## set up the minimumvariance (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
## meanvariance 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 nonBayesian monomvn, see the
## examples in the bmonomvn help file

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.