[R] fitting distributions with R

[Thomas Lumley]

On Tue, 6 Sep 2005, Ted.Harding at nessie.mcc.ac.uk wrote:

>
> However, a related method available for 'optim' is "L-BFGS-B"
> "which allows _box constraints_, that is each variable can be
> given a lower and/or upper bound. The initial value must satisfy
> the constraints." This can be set in a parameter for 'mle'.

These box constraints are really designed for situations where the 
boundary is a valid parameter value (so you are really doing constrained 
estimation) rather than situations where the boundary is an artifact of 
parameterisation.

> This interesting saga, provoked by Nadja's query, now raises
> an important general question: Given the simplicity of the
> problem, why is the use of 'mle' so unexpectedly problematic?
>

The problem is simple only in that it is one-dimensional, and optim() 
doesn't take advantage of this.  It is poorly scaled: since the starting 
value is 0.1, the maximum is at 0.00006, and there is a singularity at 0, 
it would be helpful to specify the parscale control option to optim.

The other problem is that we are using finite-difference approximations to 
the derivatives. These are bound to perform badly near the singularity at 
zero, especially in a badly scaled problem.  There is a bug in that 
L-BFGS-B doesn't respect the bounds in computing finite-differences, but 
this is not going to be easy to fix (there was recent discussion on 
r-devel about this).

If I remove the singularity by defining

> lll
function(beta) if(beta<0) 1e6 else ll(beta)

and specify parscale, I get
> est

Call:
mle(minuslogl = lll, start = list(beta = 0.01), control = list(parscale = 
1e-05))

Coefficients:
         beta
6.767725e-05

(Any parscale below 0.01 will give basically the same answer).


Incidentally, the trace output may look as if it is oscillating, but that 
is partly an artifact of the line search that BFGS uses.  The last few 
printed loglikelihoods are
[1] 254.4226
[1] 254.4226
[1] 543.2361
[1] 542.5717


Finally, as I noted earlier, this isn't really a constrained estimation 
problem, it is a problem of a function defined on an open interval with a 
singularity at one end.  In this case (in contrast to real constrained 
estimation problems) it might well be sensible to reparametrize.  mle() 
then works with no problems.

 	-thomas