I would like to fit a complex nonlinear regression model: basically, I have a complex computer code which has an input vector $\mathbf{x}$, a vector of calibration parameters $\boldsymbol{\theta}$ and produces an output $s=f(\mathbf{x},\boldsymbol{\theta})$. I also have iid observations $y$ for this output, which I model as contaminated by Gaussian, zero mean noise:
$$y = f(\mathbf{x},\boldsymbol{\theta})+\epsilon$$
Let $D=\{\mathbf{x}_i,y_i\}_{i=1}^N$ be the data set of my observations. I would like to find the vector $\boldsymbol{\theta}$ which maximizes the likelihood of observing $D$. Because of the assumption I made on the noise, maximizing the likelihood is equivalent to minimizing the residual sum of squares. This means which I can find $\boldsymbol{\theta}$ by means of Nonlinear Least Squares.
I need to use R for this project. In R, nls
is available to implement NLS, which I think implements the Newton algorithm, but I've had convergence issues using it. Also, the interface is (for me) a bit awkward.
Are there other methods for Nonlinear Least Square regression, which are readily available in R, and have better convergence properties?
EDIT I already received great suggestions (nsl2
, nlstools
and nlmrt
) which may hopefully solve my issues. Anyway, I'm providing some extra details on my case. $f(\mathbf{x},\boldsymbol{\theta})$ is a black box, i.e., it's the output of a computer code for which I don't have the source. Thus I can't provide you with an analytical, explicit expression for $f(\mathbf{x},\boldsymbol{\theta})$. I can tell you that the code computes damping and stiffness for a particular kind of turbomachinery seal, by numerically integrating the 1D isothermal compressible Navier-Stokes equations in an annular geometry. The code has been tested and validated multiple times in the course of years, so it's reliable.
Parametrization: to respond to @Roland's comment, the code is an executable with an input file, which contains many input values. Among these values, with the support of seal experts I chose some values to be predictors (vector $\mathbf{x}$) and some others to be calibration parameters (vector $\boldsymbol{\theta}$). Of course, since I don't have access to the code source, I cannot be sure that my parameters will be identifiable. But I think there must be a way to at least get an idea: I cannot open the code, but I can run it multiple times, for arbitrary values of $\mathbf{x}$ and $\boldsymbol{\theta}$ (as long as the code converges :) Can this knowledge be used to address @Roland's points? I.e.,
- whether the data support the problem (I don't understand this point)
- if I need to reparametrize (this is clear and I think it's related to identifiability)
- if I have a partial linear problem (what's this?)
EDIT2 sadly I don't have access to gradients. It's a '70 code, and writing adjoint codes was not at all common in the engineering community at that time.
nls
doesn't converge and the reason is not bad starting values, you should usually spend more time considering the maths behind your problem. Maybe your data doesn't support the model, maybe you need to re-parameterize the model, maybe you have a partial linear problem and can use an algorithm utilizing that, maybe ... $\endgroup$nls2
among my search terms on at least one of my searches though, since I recall thinking that would find some of the other packages as well. I may have searched on "starting" or "convergence"; at least those were what I had planned to search on but I don't recall how far I got through my list of things to search before I found enough to talk about. ... on the novelty of your problem, I'm not sure how you even specify the model to nls in that case. $\endgroup$nls
and if that works you are done. If that gives you problems try the approach in the following link: stackoverflow.com/questions/42511278/… $\endgroup$