I have an experimental data set $D=\{\mathbf{x}_i,y_i\}_{i=1}^N$, and a computer code with inputs $\mathbf{x}$ and calibration parameters $\boldsymbol{\theta}$, returning a value $s=f(\mathbf{x},\boldsymbol{\theta})$. Assuming Gaussian iid errors, I want to calibrate my code using the data set $D$ using nls
or similar, more advanced functions in R. If you need more details, you can find them in this question. Note that since my model is implicitly defined by my computer code and not by an explicit formula
, any nonlinear least squares function must use numerical, not analytical, Jacobians.
Now, the problem is that for some vector $\boldsymbol{\theta}^{(j)}$, the computer code doesn't converge at all points of $D$. It converges for some data points, but it doesn't converge for others. From a theoretical point of view, I believe the likelihood $\mathcal{L}(D;\boldsymbol{\theta}^{(j)})$ is not even defined.
Because of various issues, I finally switched to using nonlinear least squares functions such asnlfb
from package nlmrt
and nls.lm
from package minpack.lm
, which do not require a formula
argument, but take a function resfn
which takes in input the parameter vector $\boldsymbol{\theta}$ and return the vector of residuals $\mathbf{r}=(y_1-f(\mathbf{x}_1,\boldsymbol{\theta}),\dots,y_N-f(\mathbf{x}_N,\boldsymbol{\theta}))$. Suppose that for a certain value $\boldsymbol{\theta}^{(j)}$ proposed by the NLS algorithm, my computer code is not converging at points $\mathbf{x}_{i1},\dots,\mathbf{x}_{im}$. Now, what should I return as the corresponding components $r_{i1},\dots,r_{im}$ of the residual vector $\mathbf{r}$?
- Should I return
NA
? I don't think that would be handled bynlfb
or bynls.lm
. - Should I return an arbitrary, very large number $z$? I think this would "encourage" the NLS algorithm to move away from similar parameter values, and move towards other regions of the parameter space, where hopefully my code would converge on all of $D$. Sounds very ad-hoc.
I tried dealing with this issue empirically: I removed those few points in $D$, which correspond to conditions which are not very interesting for calibration and which are also known to pose convergence issues to my computer code. Since I have $N$=319 experimental points, and only 3 parameters to calibrate (for now!), deleting a few points wasn't an issue. However, the NLS regression is still not converging, and I need to solve this issue. Can you give me some tips?
nlfb
is Nash variant of Marquardt nonlinear least squares solution via qr linear solver. Thenls.lm
function provides an R interface to lmder and lmdif from the MINPACK library, for solving nonlinear least-squares problems by a modification of the Levenberg-Marquardt algorithm, with support for lower and upper parameter bounds. The implementation can be used via nls-like calls using the nlsLM function. Well, yuk. Maybe think better methods for a start. $\endgroup$