R devel - Aug 2008 - Numerical optimisation and "non-feasible" regions

Dear list,

I'm currently writing a C code to compute the (composite) likelihood - 
well this is done but not really robust. The C code is wrapped in an R 
one which call the optimizer routine - optim or nlm. However, the 
fitting procedure is far from being robust as the parameter space 
depends on the parameter - I have a covariance matrix that should be a 
valid one for example.

Currently, I set in my header file something like #define MINF -1.0e120 
and test if we are  in a non-feasible region, then setting the 
log-composite likelihood to MINF. The problem I see with this approach 
is that for a quite large non-feasible region, we have a kind of plateau 
where the  log-composite likelihood  is constant and may have potential 
issues with the optimizer. The other issue is that the gradient is now 
badly estimated using finite-differences.

Consequently, I'm not sure this is the most relevant approach as it 
seems that (especially the BFGS method, probably due to the estimation 
of the gradient) the optimization is really sensitive to this
"strategy"
and fails (quite often).

As I'm (really) not an expert in optimization problems, do you know good 
ways to deal with non-feasible regions? Or do I need to reparametrize my 
model so that all parameters belong to $\mathbb{R}$ - which should be 
not so easy...

Thanks for your expertise!
Best,
Mathieu

-- 
Institute of Mathematics
Ecole Polytechnique F?d?rale de Lausanne
STAT-IMA-FSB-EPFL, Station 8
CH-1015 Lausanne   Switzerland
http://stat.epfl.ch/
Tel: + 41 (0)21 693 7907

Mathieu Ribatet <mathieu.ribatet <at> epfl.ch> writes:
> 
> Dear list,
> 
> I'm currently writing a C code to compute the (composite) likelihood - 
> well this is done but not really robust. The C code is wrapped in an R 
> one which call the optimizer routine - optim or nlm. However, the 
> fitting procedure is far from being robust as the parameter space 
> depends on the parameter - I have a covariance matrix that should be a 
> valid one for example.
  
  One reasonably straightforward hack to deal with this is
to add a penalty that is (e.g.) a quadratic function of the
distance from the feasible region, if that is reasonably
straightforward to compute -- that way your function will
get gently pushed back toward the feasible region.

  Ben Bolker