# Linear regression best polynomial (or better approach to use)?

Any ideas on other polynomials I could successfully use for applying regression? My goal is a solution that has fit error strictly based on the noise. Is this possible since it is a bell-like curve? The tails can be long, short or non-existent. Am I looking for the impossible? Please note it has been a while since I have plunged this far into linear regression.

Anyway, my data points do not fit the $ax+bx^2+c$ polynomial well enough. I would like a replacement for this polynomial to apply regression (or different approach). In the following example, each point misses the fit curve by $20\%$ of the $Y$ range on average (of the data points). I looked at splines, but I do not know how to apply the regression since they look piece-wise.

Example $Y$ values (evenly spaced in $X$). It looks kind of like a bell curve. I would be applying 5 to 50 data points to determine the polynomial coefficients. Ultimately, I am interested in the best $X$ location for the peak based on the data points.

 X      Y
-40    -21142.1111111111
-30    -21330.1111111111
-20    -12036.1111111111
-10      7255.3888888889
0     32474.8888888889
10     32474.8888888889
20      9060.8888888889
30    -11628.1111111111
40    -15129.6111111111

• What are those numbers at the end of your question supposed to mean? – Maarten Buis Sep 16 '13 at 16:47
• Can you provide some context to clarify your question? I'm not sure what you are asking about. Eg, you ask about "other polynomials" before introducing the referent polynomial that you don't like. – gung - Reinstate Monica Sep 16 '13 at 16:53
• As the others, I am confused. What data have you got? What is the DV and what are the IVs? – Peter Flom Sep 16 '13 at 17:19
• Thanks for the questions and reformatting. The values at the end are Y values (spaced in X). They are put into the regression to determine the coefficients for: ax + bx^2 + c. Actual X values are -40, -30, -20, -10, 0, 10, 20, 30, 40 respectively (i.e. in order Y values appear). The X values are telescope focuser positions. Peter - I am not sure about DV and IV. If this is not sufficient information, I will look at DV and IV and get back to you. Thanks to everyone! – Buck Sep 16 '13 at 17:55
• Additionally, this is used to determine the optimal focuser position for a telescope system (i.e. using an electronic focuser and camera). Y values are brightness from real images taken at each of the focuser positions (i.e. at each X value). The goal is to find the optimal focuser position which maximizes brightness (e.g. the X location at the peak). – Buck Sep 16 '13 at 17:57

Because brightness is a response with independent random error and it is expected to taper off with distance from the optimal point according to a Gaussian function, a quick nonlinear regression ought to do a good job.

The model is

$$y = b + a \exp\left(-\frac{1}{2}\left(\frac{x-m}{s}\right)^2\right) + \varepsilon$$

where $\varepsilon$ represents the errors in measuring the brightness, modeled here as random quantities. The peak occurs at $m$; $s\gt 0$ quantifies the rate at which the curve tapers off; $a\gt 0$ reflects the overall magnitudes of the relative $y$ values, and $b$ is a baseline.

Let's try it with the sample data (using R). By including the middle ($m$) among the parameters, the software will automatically output its estimate and a standard error for it:

y <- c(-190279, -191971, -108325, 65298, 292274, 292274, 81548, -104653, -136166)/9
x <- (-4:4)*10
#
# Define a Gaussian function (of four parameters).
f <- function(x, theta)  {
m <- theta[1]; s <- theta[2]; a <- theta[3]; b <- theta[4];
a*exp(-0.5*((x-m)/s)^2) + b
}
#
# Estimate some starting values.
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; b.0 <- min(y); a.0 <- (max(y)-min(y))
#
# Do the fit.  (It takes no time at all.)
fit <- nls(y ~ f(x,c(m,s,a,b)), data.frame(x,y), start=list(m=m.0, s=s.0, a=a.0, b=b.0))
#
# Display the estimated location of the peak and its SE.
summary(fit)$parameters["m", 1:2]  The tricky part with nonlinear fits usually is finding good starting values for the parameters: this code shows one (crude) approach. Its output,  Estimate Std. Error 5.3161940 0.4303487  gives the estimate of the peak ($5.32$) and its standard error ($0.43$). It's always a good idea to plot the fit and compare it to the data: par(mfrow=c(1,1)) plot(c(x,0),c(y,f(coef(fit)["m"],coef(fit))), main="Data", type="n", xlab="x", ylab="Brightness") curve(f(x, coef(fit)), add=TRUE, col="Red", lwd=2) points(x,y, pch=19)  That's what we expected: the data appear to fit the Gaussian pretty well. For a more incisive look at the fit, plot the residuals: plot(x, resid(fit), main="Residuals")  You want to check that most residuals are as small as the (known?) variation in the brightness measurement and that there is no important trend or pattern in them. We might be a little concerned about the high residual at$x=40$, but rerunning the procedure with this last data point removed does not appreciably change the estimate of$m$(which is now$5.25$with a standard error of$0.17$, not distinguishable from the previous estimate). The new residuals bounce up and down, tend to get smaller as$x$gets larger, but otherwise tend to be less than$1000$or so in absolute value: there's no sign here that more effort is needed to pin down$m$. • Perfect as far as what I was hoping for. It is amazing how well the data matches the curve. Thank you! – Buck Sep 16 '13 at 21:31 • @whuber You wrote$a>0$and$s>0$, does the R code that you wrote respect these constraints? – Alessandro Jacopson Oct 8 '14 at 8:36 • @uvts That's a good question.$s\gt 0$is not a constraint because the model for$-s$is the same as that for$s$: this is just a way of assuring identifiability of$s$. Rather than constrain$a$to be non-negative (which could be done) I find it both more informative and more general to allow the optimizer to fit the best unconstrained value of$a$it can. – whuber Oct 8 '14 at 14:31 I assume the column of values you gave correspond to some kind of time series, and there is an implicit "time" column with evenly spaced values that you've not mentioned. Given that, the question remains if you intend to fit a polynomial curve for prediction or forecasting. If the former, then the use of splines can be achieved using adaptive bases and cross validation to determine optimal breakpoints. Splines estimate piecewise polynomial trends, depending on the polynomial degree you specify (and the number of breakpoints). On the other hand, forecasting does not warrant the use of splines with breakpoints because the trend cannot be extrapolated beyond where the next breakpoint would lie had you observed the data. In either case, are you determining the$R^2$value using external validation with an independent dataset? If so, then how did you arrive at a choice of 80% variation to determine an adequate predictive model? It seems arbitrary to me and you're more likely to fit noise with little generalization by doing so. • Yes, this is a time series for prediction. The determined coefficients are a=200.94, b=-27.56, c=9189.85. My measure of fit is a percent of the range for the data points specified. In this case, it is 20% and is computed by traversing each data point and getting the absolute difference to the determined curve. These are averaged and reported as a percent of the data point range (i.e. max - min of all data points). Feel free to poke holes at this. I am a developer trying to get some good results. I expect R^2 may be better but this is pretty easy to visualize. – Buck Sep 16 '13 at 18:37 • Is it so hard to cut off the other 12 decimals? Pearson$R^2\$ is already standardized to the variance of the data, so no need to point this out. You should be using the variance of residuals to determine how good this fit is. – AdamO Sep 16 '13 at 18:39
• Oh yes, sorry if you saw any previous partial comment. I was pressing enter which submitted it early. I will get use to this. Sorry about the extra digits - now edited. Thanks. – Buck Sep 16 '13 at 18:39
• Got it. Can I use splines if I do not know which of my data points match up with which pieces (since piecewise). If so, I will the research. Thanks! – Buck Sep 16 '13 at 18:58
• Yes, you can use cross validation with a number of possible spline breakpoints to determine where the best cut off lies empirically. However, beware that this can affect generalizability and hurt your prediction in other datasets. – AdamO Sep 16 '13 at 20:31

If your data should be a bell curve, you should fit that instead of the quadratic you have in the question. However, you can play with "other polynomials", for example in Python using mean squared error minimization with numpy.polyfit:

import numpy as np
import pylab as pl
# generate bell curve data
X = np.sort((10 * np.random.rand(50, 1)-5), axis=0)
y = ((np.exp(-(X**2)/2))/(2*np.pi)).ravel()
y[::5] += (0.1 * np.random.rand(len(y)/5))
x = X.ravel()
# do the fit with 4 inflections and get the coefficients
z = np.polyfit(x,y,4)
# get a callable
p = np.poly1d(z)
# plot
hf = pl.figure()
xlab = xrange(0,len(x))
ax.scatter(x, y, c='k', label='data')
ax.plot(X,p(X), c='g', label='fit')
ax.set_xticklabels(xlab)
ax.set_xlabel('time')
ax.set_ylabel('brightness')
ax.set_title('Polyfit')
ax.legend()
pl.show()


Output:

• You really don't want to fit polynomials to a Gaussian: the tails, which provide next to no information about the location of the peak, have too much influence in the fitting. – whuber Sep 16 '13 at 20:30