How to curve-fit the lower envelope of random sequence?

Question

I'm more or less familiar with procedures and methods to fit a curve to experimental data, and I have done this many times using Matlab. However this time I have a problem that I'm not sure how to solve and I'm stuck.

I have some experimental data which I have plotted in the figure below (blue curve).

The data seem to be lowerbounded by some kind of smooth logarithmic-like curve like the red line shown in the figure (I have plotted this red curve manually using Paint just to illustrate what I mean).

I have several mathematical expressions that could be good candidates to fit the red line, but I'm not sure how to fit them in the way shown in the figure. Standard data fitting procedures use the whole data set and provide a curve somewhere in the middle of the blue curve, but I need something like the red curve.

I work with Matlab and I have tried to use the envelope function (lower envelope), the absolute value of the Hilbert transform and the convex hull (convhull) but nothing works.

Any ideas or suggestions'

Answer 1

Well, let our data be $(x_i, y_i)_{i=1}^{n}$ and suppose you have a parametric model $f(x|\beta)$ with parameters $\beta$, that you want to fit to your data. A simple formulation that will help give you what you want is:

$$\begin{align*} &\min_{\beta, \tau} \sum_{i=1}^n (y_i - f(x_i|\beta)) + \gamma\sum_{j=1}^{m} \tau_j\\ &\text{subject to}\\ &\hspace{2cm} \forall i, y_i \geq f(x_i|\beta)\\ &\hspace{2cm} \forall j, \beta_j \leq \tau_j \\ &\hspace{2cm} \forall j, \beta_j \geq -\tau_j \\ &\hspace{2cm} \forall j, \tau_j \geq 0 \end{align*}$$

where ultimately $\beta$ are the model parameters you want to find, $\gamma$ is some regularization parameter which you can tune via k-fold validation, and $\tau_j$ are just added variables to enforce an $L_1$ regularization on $\beta$, if you decide you want that. The reason I chose this particular form is because if your model happens to have the form $$f(x|\beta) = \sum_{j=1}^{m} \beta_j g_j(x)$$ ie your function is linear in your model parameters, then this implies the above formulation is a linear program and can be solved efficiently.

If your function is nonlinear in $\beta$ or you want the objective to involve minimizing squared loss or using $L_2$ regularization, then the computational becomes more hairy.

Answer 2

Assuming you have a model $M(x)$ that you can evaluate at $x$ and compare to your data at each point $Y(x)$.

Try minimizing $\sum D^2$ where $D = Y(x) – M(x)$ if $Y(x) > M(x)$ and $D = P$ if $Y(x) < M(x)$ and $P$ is some imposed penalty if $Y(x) < M(x)$.

Try smaller and smaller $P$ until too many of the $Y(x)$ are below $M(x)$.