因为亮度是具有独立随机误差的响应，并且根据高斯函数预计会随着与最佳点的距离而逐渐减小，所以快速非线性回归应该做得很好。

模型是

在哪里$\varepsilon$表示测量亮度的误差，这里建模为随机量。峰值出现在$m$;$s\gt 0$量化曲线逐渐变细的速率；$a\gt 0$反映相对的整体大小$y$价值观，以及$b$是基线。

让我们用示例数据（使用R）来尝试一下。通过包括中间 ($m$) 在这些参数中，软件会自动输出它的估计值和一个标准误差：

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]

非线性拟合的棘手部分通常是为参数找到好的起始值：此代码显示了一种（粗略的）方法。它的输出，

  Estimate Std. Error 
 5.3161940  0.4303487 

给出峰值的估计值 ($5.32$) 及其标准误 ($0.43$）。绘制拟合并将其与数据进行比较总是一个好主意：

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)

这正是我们所期望的：数据似乎非常适合高斯分布。为了更深入地了解拟合，绘制残差：

plot(x, resid(fit), main="Residuals")

您想检查大多数残差是否与亮度测量中的（已知？）变化一样小，并且其中没有重要的趋势或模式。我们可能有点担心高残差$x=40$，但在删除最后一个数据点的情况下重新运行该过程不会明显改变$m$（现在是$5.25$标准误为$0.17$，与之前的估计没有区别）。新的残差上下反弹，随着$x$变大，但在其他方面往往小于$1000$绝对值左右：这里没有迹象表明需要更多的努力来确定$m$.

我假设您提供的值列对应于某种时间序列，并且有一个隐含的“时间”列，其中包含您没有提到的均匀间隔的值。

鉴于此，如果您打算拟合多项式曲线进行预测或预测，问题仍然存在。如果是前者，那么可以使用自适应基础和交叉验证来确定最佳断点，从而实现样条的使用。样条曲线估计分段多项式趋势，具体取决于您指定的多项式次数（和断点数）。

另一方面，预测并不保证使用带有断点的样条曲线，因为如果您观察到数据，则无法将趋势外推到下一个断点所在的位置之外。

无论哪种情况，您是否确定$R^2$使用具有独立数据集的外部验证的价值？如果是这样，那么您是如何选择 80% 的变异来确定适当的预测模型的？在我看来，这似乎是任意的，并且您更有可能通过这样做来适应噪音而几乎没有泛化。

如果您的数据应该是钟形曲线，您应该拟合它而不是问题中的二次曲线。但是，您可以使用“其他多项式”，例如在 Python 中使用 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()
# add some noise
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()
ax = hf.add_subplot(1,1,1)
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()

输出：