[R] question about nlminb
Spencer Graves
spencer.graves at pdf.com
Sun Apr 13 22:40:37 CEST 2008
It's very well known that if a random vector X has a finite mean
mu and covariance Sig, and Y = A X, then
(1) EY = A %*% mu
and
(2) cov(Y) = A %*% Sig %*% t(X)
= tcrossprod(A %*% Sig, A)
Expression (1) says that mathematical expectation is a linear
operator. Expression (2) follows by application of (1). See any good
reference that discusses expected values. I didn't find a web site with
a 2 minute search, but I would expect to find it in C. R. Rao (1973)
Linear Statistical Inference and Its Application, 2nd ed. (Wiley) and in
many other places.
Best Wishes,
Spencer
John Pitchard wrote:
> 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.
>
More information about the R-help
mailing list