[R] fitting distributions with R

Thomas Lumley tlumley at u.washington.edu
Tue Sep 6 18:50:48 CEST 2005

```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

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

```