[R] help with the maxBHHH routine
Andrew Robinson
A.Robinson at ms.unimelb.edu.au
Thu May 5 01:28:41 CEST 2011
Hi Rohit,
actually, the request for simple reproducible code means that you have
to find the simplest possible representation of the problem.
What happens if you simplify the observation level gradient and the
likelihood function? Eg to trivial examples? If you still get the
error, then simplify it futher. If you get the error with the
simplest possible problem, then share it. If you don't , then try to
figure out what the changes were that resolved the problem, and scale
those back up to your original problem.
Does that make sense?
Cheers
Andrew
On Thu, May 05, 2011 at 03:22:55AM +0530, Rohit Pandey wrote:
> Hi Andrew, Ravi and Arne,
>
> Thank you so much for your prompt replies. I see that all of you mention
> the need for simple, reproducible code. I had thought of doing this, but
> the functions I was using for the observation level gradient and
> likelihood function were very long. I will paste them below here.
>
> Also, sorry for the ambiguity with the "1000's of observations and 821
> parameters" on the one hand and the 10 * 2 matrix on the other. The latter
> is a toy data set and the former is the real data set I ultimately hope to
> apply this routine to once it works. Also, sorry for not mentioning the
> fact that the maxBHHH function I am using is from the maxLik package
> (thanks, Ravi for pointing out).
> So, the code that is giving me the errors is:
> maxBHHH(logLikALS4,grad=nuGradientC4,finalHessian="BHHH",start=prm,iterlim=2)
> and
> maxBHHH(logLikALS4,grad=nuGradientC4,finalHessian="BHHH",start=prm,iterlim=2)
> Where nuGradientC4 returns a 2*10 matrix and nuGradientC5 a 10*2 matrix
> (there are 10 parameters and 2 observations).
> I have attached the required functions in the .R file.
> These make for some pretty long code, but all you have to do is either
> load the file or paste the contents into your R console (and maybe see
> that they're returning what they're supposed to). I'm sorry I couldn't
> think of a way to come up with a shorter version of this code (I tried my
> best).
>
> Once you load the file, you should see the following:
>
> #The observation level likelihood function
> > logLikALS4(prm)
> 1 2
> -0.6931472 -0.6931472
>
> #The observation level gradients
> > nuGradientC4(prm)
> 1 2 3 4 5 6 7 8 9 10
> 2 -0.3518519 0.3518519 0.0000000 0 -0.1481481 -0.1666667 0.1481481
> 0.1666667 0.0000000 0.0000000
> 4 0.0000000 -0.3518519 0.3518519 0 0.0000000 0.0000000 -0.1666667
> -0.1481481 0.1666667 0.1481481
> Warning messages:
> 1: In [1]is.na(x) : [2]is.na() applied to non-(list or vector) of type
> 'NULL'
> 2: In [3]is.na(x) : [4]is.na() applied to non-(list or vector) of type
> 'NULL'
>
> > nuGradientC5(prm)
> 2 4
> 1 -0.3518519 0.0000000
> 2 0.3518519 -0.3518519
> 3 0.0000000 0.3518519
> 4 0.0000000 0.0000000
> 5 -0.1481481 0.0000000
> 6 -0.1666667 0.0000000
> 7 0.1481481 -0.1666667
> 8 0.1666667 -0.1481481
> 9 0.0000000 0.1666667
> 10 0.0000000 0.1481481
> Warning messages:
> 1: In [5]is.na(x) : [6]is.na() applied to non-(list or vector) of type
> 'NULL'
> 2: In [7]is.na(x) : [8]is.na() applied to non-(list or vector) of type
> 'NULL'
>
> Ignore the warning messages.
>
> The errors are:
>
> >
> maxBHHH(logLikALS4,grad=nuGradientC4,finalHessian="BHHH",start=prm,iterlim=2)
> Error in checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
> :
> the matrix returned by the gradient function (argument 'grad') must have
> at least as many rows as the number of parameters (10), where each row
> must correspond to the gradients of the log-likelihood function of an
> individual (independent) observation:
> currently, there are (is) 10 parameter(s) but the gradient matrix has only
> 2 row(s)
> In addition: Warning messages:
> 1: In [9]is.na(x) : [10]is.na() applied to non-(list or vector) of type
> 'NULL'
> 2: In [11]is.na(x) : [12]is.na() applied to non-(list or vector) of type
> 'NULL'
>
> and:
>
> >
> maxBHHH(logLikALS4,grad=nuGradientC5,finalHessian="BHHH",start=prm,iterlim=2)
> Error in gr[, fixed] <- NA : (subscript) logical subscript too long
> In addition: Warning messages:
> 1: In [13]is.na(x) : [14]is.na() applied to non-(list or vector) of type
> 'NULL'
> 2: In [15]is.na(x) : [16]is.na() applied to non-(list or vector) of type
> 'NULL'
>
> Again, thanks for your patience and help.
>
> Rohit
> On Wed, May 4, 2011 at 4:44 AM, Andrew Robinson
> <[17]A.Robinson at ms.unimelb.edu.au> wrote:
>
> I suggest that you provide some commented, minimal, self-contained,
> reproducible code.
>
> Cheers
>
> Andrew
> On Wed, May 04, 2011 at 02:23:29AM +0530, Rohit Pandey wrote:
> > Hello R community,
> >
> > I have been using R's inbuilt maximum likelihood functions, for the
> > different methods (NR, BFGS, etc).
> >
> > I have figured out how to use all of them except the maxBHHH function.
> This
> > one is different from the others as it requires an observation level
> > gradient.
> >
> > I am using the following syntax:
> >
> >
> maxBHHH(logLik,grad=nuGradient,finalHessian="BHHH",start=prm,iterlim=2)
> >
> > where logLik is the likelihood function and returns a vector of
> observation
> > level likelihoods and nuGradient is a function that returns a matrix
> with
> > each row corresponding to a single observation and the columns
> corresponding
> > to the gradient values for each parameter (as is mentioned in the
> online
> > help).
> >
> > however, this gives me the following error:
> >
> > *Error in checkBhhhGrad(g = gr, theta = theta, analytic =
> (!is.null(attr(f,
> > :
> > the matrix returned by the gradient function (argument 'grad') must
> have
> > at least as many rows as the number of parameters (10), where each row
> must
> > correspond to the gradients of the log-likelihood function of an
> individual
> > (independent) observation:
> > currently, there are (is) 10 parameter(s) but the gradient matrix has
> only
> > 2 row(s)
> > *
> > It seems it is expecting as many rows as there are parameters. So, I
> changed
> > my likelihood function so that it would return the transpose of the
> earlier
> > matrix (hence returning a matrix with rows equaling parameters and
> columns,
> > observations).
> >
> > However, when I run the function again, I still get an error:
> > *Error in gr[, fixed] <- NA : (subscript) logical subscript too long*
> >
> > I have verified that my gradient function, when summed across
> observations
> > gives the same results as the in built numerical gradient (to the 11th
> > decimal place - after that, they differ since R's function is
> numerical).
> >
> > I am trying to run a very large estimation (1000's of observations and
> 821
> > parameters) and all of the other methods are taking way too much time
> > (days). This method is our last hope and so, any help will be greatly
> > appreciated.
> >
> > --
> > Thanks in advance,
> > Rohit
> > Mob: 91 9819926213
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > [18]R-help at r-project.org mailing list
> > [19]https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> [20]http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> --
> Andrew Robinson
> Program Manager, ACERA
> Department of Mathematics and Statistics Tel: +61-3-8344-6410
> University of Melbourne, VIC 3010 Australia (prefer email)
> [21]http://www.ms.unimelb.edu.au/~andrewpr Fax: +61-3-8344-4599
> [22]http://www.acera.unimelb.edu.au/
>
> Forest Analytics with R (Springer, 2011)
> [23]http://www.ms.unimelb.edu.au/FAwR/
> Introduction to Scientific Programming and Simulation using R (CRC,
> 2009):
> [24]http://www.ms.unimelb.edu.au/spuRs/
>
> --
> Thanks,
> Rohit
> Mob: 91 9819926213
>
> References
>
> Visible links
> 1. http://is.na/
> 2. http://is.na/
> 3. http://is.na/
> 4. http://is.na/
> 5. http://is.na/
> 6. http://is.na/
> 7. http://is.na/
> 8. http://is.na/
> 9. http://is.na/
> 10. http://is.na/
> 11. http://is.na/
> 12. http://is.na/
> 13. http://is.na/
> 14. http://is.na/
> 15. http://is.na/
> 16. http://is.na/
> 17. mailto:A.Robinson at ms.unimelb.edu.au
> 18. mailto:R-help at r-project.org
> 19. https://stat.ethz.ch/mailman/listinfo/r-help
> 20. http://www.r-project.org/posting-guide.html
> 21. http://www.ms.unimelb.edu.au/%7Eandrewpr
> 22. http://www.acera.unimelb.edu.au/
> 23. http://www.ms.unimelb.edu.au/FAwR/
> 24. http://www.ms.unimelb.edu.au/spuRs/
--
Andrew Robinson
Program Manager, ACERA
Department of Mathematics and Statistics Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia (prefer email)
http://www.ms.unimelb.edu.au/~andrewpr Fax: +61-3-8344-4599
http://www.acera.unimelb.edu.au/
Forest Analytics with R (Springer, 2011)
http://www.ms.unimelb.edu.au/FAwR/
Introduction to Scientific Programming and Simulation using R (CRC, 2009):
http://www.ms.unimelb.edu.au/spuRs/
More information about the R-help
mailing list