Re: [R] question about nlminb

From: John Pitchard <johnpitchard_at_googlemail.com>
Date: Sun, 13 Apr 2008 21:03:39 +0100

Hi Spencer,

Thanks for your email.

Do you have a reference for generating the variance-covariance matrix from the restricted variance-covariance? Is this a well known technique?

Regards,
John

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.

list of date sections of archive