# 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);
tmp(6,5) = 1.2/(sig*sqrt(2*pi));
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.

1. 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.
2. 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 "binary" data, but rather a uniform sampled array with amplitude.
3. 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!

• Thank for the answers. I think that part of what confuses me is that simple explanations of regression smoothing and tracking tends to focus on "scalar" sensory input. You might have raw estimates of position, speed, acceleration etc in 3 dimensions + 3 rotations. But each sensor reports a scalar value (typically as a function of time). What I want to wrap my head around here is that there is no "y_hat=f(x)" sensory input. There is only a 2-d array of (noisy) "probabilities", and rather than preprocessing to produce an imperfect y_hat, I'd like to work directly on the input array. Jun 26 '20 at 13:30

Just a couple of thoughts:

• That sounds like the classical application for a Kalman filter: you have one gaussian-overlaid signal that "develops" 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 &= \arg\max\limits_\mu f_{\mu|\mathbf C}(\mu)\\ &= \arg\max\limits_\mu\;\log\left(f_{\mu|\mathbf C}\right)\\ &= \arg\max\limits_\mu\;\log \left(\prod\limits_{\mathbf c \in \mathbf C}f_{\mu|\mathbf c} \right)\\ &= \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.

You need to think 3D. You are looking at the domain and the range is the intensity.

The first step is to fix the outliers as you have them. For each cell in the plane, calculate the average of the neighbors (before their adjustments). If the cell value exceeds a specified threshold, limit the value to that threshold.

Next do four smoothing passes: Horizontal, vertical, diagonal , diagonal /. For each direction, for each strip, smooth with a symmetric filter so no phase shifts, find the peak with parabolic estimation, and create a point (x,y) where the peak occurs.

For conic sections, use this formula (I go cubic).

$$f(x,y) = Ax^2 + Bxy + Cy^2 + Dx + Ey$$

Find the best fit (A,B,C,D,E) for $$f=0$$. For a given level curve $$f(x,y)= L$$, you have a conic section in (x,y).

Obviously, if the data allows, you can also curve fit to:

$$y = Ax^4 + Bx^3 + Cx^2 + Dx + E$$