Fitting a polynomial to a ridge in an array I have an array of numerical data. Say,
%% setup an array of gaussians where the mean is quadratically swept
p = [0.15,-0.3,1];
mu = polyval(p, 1:8);
sig = 2;
tmp = (1/(sig*sqrt(2*pi)))*exp(-(1/2)*(((1:8)'-mu)/sig).^2);
%% add interference
tmp(6,5) = 1.2/(sig*sqrt(2*pi));
%% add noise
tmp = tmp + 1.5e-2*randn(size(tmp));


Assuming that the gaussians are oversampled (sample maxima are close enough to true maxima).
I would like to estimate p from the output array, p_hat.
I have been looking at it for some time, and I cannot twist my head around how to relate this to stuff that I know.

The most straight forward solution would be to take the maximum of each column, do a 1-d fitting using ordinary least squares, solve the outlier somehow. But then the column containing the outlier is not counted at all. I was hoping to be smarter than that, to improve the fit by depending on the secondary (true) peak of the interfered column.
On the face of it, this sounds related to ordinary least squares. But there are 2 independent variables (x-dim and y-dim). And the data is not a sparse point-cloud of &quot;binary&quot; data, but rather a uniform sampled array with amplitude.
The thought to traversing the array left to right (making assumptions of continuity and start or end-point) in a linear-programming kind of way occured to me, but then I would be fitting a continous line with constraints, rather than a parametric polynoma.

I prefer to not use toolboxes. Thanks for any ideas!
Just a couple of thoughts:

That sounds like the classical application for a Kalman filter: you have one gaussian-overlaid signal that &quot;develops&quot; over the columns
If you have a signal model for your signal of interest, a parametric estimator might indeed be the best choice here.
You could write a likelihood function for each observation column $\mathbf c$, like $f_{\mu|\mathbf c}(\mu)=\ldots$ that describes how likely a given $\mu$ is given the observation in each column.
Then, you could state $\hat \mu = \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)$.
If you can assume your columns to be independent: $f_{\mu|\mathbf C}=\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c}$, so that
\begin{align}\hat \mu &amp;= \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)\\
&amp;= \arg\max\limits_\mu\;\log\left(f_{\mu|\mathbf C}\right)\\
&amp;= \arg\max\limits_\mu\;\log \left(\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c} \right)\\
&amp;= \arg\max\limits_\mu\;\sum\limits_{\mathbf c \in \mathbf C}\log\;f_{\mu|\mathbf c}\text,\end{align}(the logarithm might be very handy, considering the exponential nature your signal of interest seems to have). From there, it's an optimization problem.