From: John Pitchard <johnpitchard_at_googlemail.com>

Date: Sun, 13 Apr 2008 21:03:39 +0100

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 13 Apr 2008 - 20:08:59 GMT

Date: Sun, 13 Apr 2008 21:03:39 +0100

On 10/04/2008, Spencer Graves <spencer.graves_at_pdf.com> wrote:

> Hi, John:

*> I just got the following error right after the attempt to use
**> 'rmvnorm'.
**> Error: could not find function "rmvnorm"
**>
**> I tried 'library(mvtnorm)', but the 'rmvnorm' in that package gave me
**> the following:
**> Error in rmvnorm(10000, mean = c(3, -20, -10, 3, 2), sd = c(0.1, 15, 4, :
**> unused argument(s) (sd = c(0.1, 15, 4, 0.15, 0.8), cov = c(1, 0.5, 0.5,
**> 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 1,
**> 0.5, 0.5, 0.5, 0.5, 0.5, 1))
**>
**> Then I got 87 hits to RSiteSearch("rmvnorm", 'fun').
**> Regarding, "How would I go about getting the entire variance-covariance
**> matrix?", this is really easy: Just define a 5 x 4 matrix A so the
**> 5-dimensional 'opt' you want is a constant plus A times the 4-dimensional
**> restricted 'opt'. [Please don't use the same name for two different things
**> like 'opt' in your function below. Half a century ago, when Fortran was
**> young, that was very smart programming. Today, it's primary effect it to
**> make it difficult for mere mortals to understand your code. Use something
**> like 'opt5' for one and 'opt4' for the other.]
**> Then
**>
**> Sig5 = A %*% vcov.nlminb(opt) %*% t(A).
**> I believe the A matrix you want is as follows:
**> A = matrix(NA, dim=c(5, 4))
**> A[1:4, 1:4] <- diag(4)
**> A[5, ] <- (-1)
**>
**> However, since your example was not self contained, I have not tested
**> this. Hope this helps. Spencer
**>
**>
**> John Pitchard wrote:
**>
**> > Hi Spencer,
**> >
**> > Sorry for not producing code as a worked example.
**> >
**> >
**> > Here's an example:
**> > ==================================
**> > # setting the seed number
**> > set.seed(0)
**> > # creating a correlation matrix
**> > corr <- diag(5)
**> > corr[lower.tri(corr)] <- 0.5
**> > corr[upper.tri(corr)] <- 0.5
**> >
**> > # Data for the minimisation
**> > mat <- rmvnorm(10000, mean=c(3, -20, -10, 3, 2), sd=c(0.1, 15, 4,
**> > 0.15, 0.8), cov=corr)
**> >
**> > obj.fun <- function(opt, mat) {
**> > opt <- c(opt, 1-sum(opt))
**> > LinearComb <- mat%*%opt
**> > obj <- -min(LinearComb)
**> > obj
**> > }
**> >
**> > opt <- nlminb(rep(0,4), lower=rep(-3, 4), upper=rep(3, 4), obj.fun,
**> mat=mat)
**> > opt.x <- opt$parameters
**> > opt.x <- c(opt.x, 1-sum(opt.x))
**> >
**> > # using vcov.nlminb from the MASS library to obtain the covariance matrix
**> > vcov.nlminb(opt)
**> > ====================================
**> >
**> >
**> > I have a variance-covariance matrix for 4 of the elements in the
**> > vector but not the last component. How would I go about getting the
**> > entire variance-covariance matrix?
**> >
**> > Thanks in advance for any help.
**> >
**> > Regards,
**> > John
**> >
**> >
**> >
**> >
**> > On 09/04/2008, Spencer Graves <spencer.graves_at_pdf.com> wrote:
**> >
**> >
**> > > Have you considered optimizing over x1 = x[1:(length(x)-1]? You
**> could feed a wrapper function 'f2(x1, ...)' that computes xFull = c(x1,
**> 1-sum(x1)) and feeds that to your 'fn'.
**> > > If this makes sense, great. Else, if my answer is not useful, be so
**> kind as to PLEASE do read the posting guide
**> http://www.R-project.org/posting-guide.html and provide
**> commented, minimal, self-contained, reproducible code.
**> > > Spencer
**> > >
**> > > John Pitchard wrote:
**> > >
**> > >
**> > >
**> > > > Dear All,
**> > > >
**> > > > I wanted to post some more details about the query I sent to s-news
**> last
**> > > > week.
**> > > >
**> > > > I have a vector with a constraint. The constraint is that the sum of
**> the
**> > > > vector must add up to 1 - but not necessarily positive, i.e.
**> > > >
**> > > > x[n] <- 1 -(x[1] + ...+x[n-1])
**> > > >
**> > > > I perform the optimisation on the vector x such that
**> > > >
**> > > > x <- c(x, 1-sum(x))
**> > > >
**> > > > In other words,
**> > > >
**> > > > fn <- function(x){
**> > > > x <- c(x, 1 - sum(x))
**> > > > # other calculations here
**> > > > }
**> > > >
**> > > > then feed this into nlminb()
**> > > >
**> > > > out <- nlminb(x, fn)
**> > > > out.x <- out$parameters
**> > > > out.x <- c(out.x, 1 - sum(out.x))
**> > > > out.x
**> > > >
**> > > > I would like to calculate standard errors for each of the components
**> of x.
**> > > > Is this possible by outputing the Hessian matrix? Furthermore, how
**> would I
**> > > > calculate this for the last component (if this is indeed possible)
**> which has
**> > > > the restriction (i.e. 1-sum(out.x))?
**> > > >
**> > > > Any help would be much appreciated.
**> > > >
**> > > > Regards,
**> > > > John
**> > > >
**> > > > [[alternative HTML version deleted]]
**> > > >
**> > > > ______________________________________________
**> > > > R-help_at_r-project.org mailing list
**> > > > https://stat.ethz.ch/mailman/listinfo/r-help
**> > > > PLEASE do read the posting guide
**> http://www.R-project.org/posting-guide.html
**> > > > and provide commented, minimal, self-contained, reproducible code.
**> > > >
**> > > >
**> > > >
**> > > >
**> > >
**> >
**> > ______________________________________________
**> > R-help_at_r-project.org mailing list
**> > https://stat.ethz.ch/mailman/listinfo/r-help
**> > PLEASE do read the posting guide
**> http://www.R-project.org/posting-guide.html
**> > and provide commented, minimal, self-contained, reproducible code.
**> >
**> >
**> >
**>
*

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 13 Apr 2008 - 20:08:59 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 13 Apr 2008 - 21:30:29 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*