[R] maxLik: different estimations from different packages

[Spencer Graves]

On 5/10/2013 5:22 AM, Ott Toomet wrote:
> Hi,
> maybe we changed the convergence criteria a little bit?  I cannot remember,
> but we have played around with those a few times.  You may try those,
> perhaps it helps a bit.
>
> Another recommendation---as you do not supply analytic gradient, you may
> want use BHHH method instead (if this applies to your problem), or another
> one that does not need Hessian.  Numeric approximation of Hessian is pretty
> error-prone..


       You can try all the different "method" options for "maxLik", pick 
the best answer, iterate once to make sure you can't get something 
better.  If you want a hessian and don't get one from the best answer, 
you can feed the best answer to the "hessian" function to get that.


       Also, if you get a singular hessian, that says that you cannot 
estimate as many parameters as you have.  If you have p parameters and 
the hessian is of rank p0 < p, that says you can only estimate p0 
parameters independently.  You can use the "fnSubset" function to fix 
(p-p0) parameters at specific values and see what you get from varying 
the others.


       Estimating hessians can be error-prone whether you use numeric or 
analytic derivatives.  The "compareDerivatives" function can help you 
check analytic derivatives for software bugs.


       Hope this helps.
       Spencer

> Best,
> Ott
>
>
> On Fri, May 10, 2013 at 2:40 PM, Arne Henningsen
> <arne.henningsen at gmail.com>wrote:
>
>> Dear Alfonso
>>
>> On 10 May 2013 12:51,  <alfonso.carfora at uniparthenope.it> wrote:
>>> we are computing maximum likelihood estimations using maxLik package and
>> we
>>> realized that the results of the estimation depend on the version of the
>>> package installed
>>>
>>> for example if we try to estimate this function with an old version of
>>> maxLik under R 2.13.1 (32 bit version installed 2 years ago):
>>>
>>>
>>> L<-function (param) {b0t<-param[1]
>>>                       p1t<-param[2]
>>>                       p2t<-param[3]
>>>                       p3t<-param[4]
>>>                       p4t<-param[5]
>>>
>>> for(i in 17:T) {n[i,]<- b0t + p1t*a[i-1] + p2t*sum(a[(i-4):(i-1)]) +
>>> p3t*(sum(a[(i-8):(i-1)])) + p4t*(sum(a[(i-16):(i-1)]))
>>>                 m[i,]<-exp(n[i])/(1+exp(n[i]))
>>>                 ll[i-16,]<-a[i]*log(m[i])+(1-a[i])*log(1-m[i]) }
>>>                       sum(ll)}
>>> b2<-maxLik(L, start=c(-2.8158,5,-1,0.3213,-0.3112))
>>>
>>>
>>> we obtain this results:
>>>
>>> summary(b2)
>>>
>>> Maximum Likelihood estimation
>>> Newton-Raphson maximisation, 16 iterations
>>> Return code 2: successive function values within tolerance limit
>>> Log-Likelihood: -38.11285
>>> 5  free parameters
>>> Estimates:
>>>        Estimate Std. error t value   Pr(> t)
>>> [1,]  -2.81578    0.43548 -6.4660 1.007e-10 ***
>>> [2,]  50.50024   13.17046  3.8344 0.0001259 ***
>>> [3,] -11.53344    3.31075 -3.4836 0.0004947 ***
>>> [4,]   0.32130    0.42978  0.7476 0.4547096
>>> [5,]  -0.31121    0.23245 -1.3388 0.1806280
>>> ---
>>> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>>> --------------------------------------------
>>>
>>> NB: a is a binary time series
>>>
>>> if we try to estimate the same function using the last version of maxLik
>>> under R.3.0 (64 bit latest version) the estimation do not converge and
>> this
>>> is the error message:
>>>
>>> Iteration 15
>>> Parameter:
>>> [1]  -2.8146429  51.3042554 -11.7373313   0.3245214  -0.3125767
>>> Gradient:
>>> [1] NaN NaN NaN NaN NaN
>>> Errore in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
>>> hessOrig = hess,  :
>>>    NA in gradient
>>>
>>> What causes this?
>> It could be that the NaNs in the gradients are caused by rounding
>> errors and/or approximation errors in the numerical
>> (finite-difference) derivatives when using R.3.0 64 bit, because
>> different hardware and different software versions (e.g. R,
>> mathematical libraries, OS) could lead to different rounding errors.
>> In this case, the specification of a function that returns analytical
>> gradients could solve the problem.
>>
>> If this does not solve the problem and you cannot find out the reason
>> for the NaNs in the analytical gradients yourself, please provide a
>> reproducible example so that we could help you with this.
>>
>> Please note that you could also ask questions regarding the maxLik
>> package via a forum at maxLik's R-Forge site:
>>
>> https://r-forge.r-project.org/projects/maxlik/
>>
>>
>> ... and please do not forget to cite maxLik in your publications :-)
>>
>> Best regards,
>> Arne
>>
>>
>> --
>> Arne Henningsen
>> http://www.arne-henningsen.name


-- 
Spencer Graves, PE, PhD
President and Chief Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567
web:  www.structuremonitoring.com