寻找一条曲线来匹配数据

Question

我正在寻找一种非线性曲线拟合程序（可能最容易在R或Python中找到，但我也可以考虑其他编程语言），该程序可以采用x、y数据并对其进行曲线拟合。

我应该能够通过一个字符串来指定我想要拟合的表达式类型。

示例：

"A+B*x+C*x*x"
"(A+B*x+C*x*x)/(D*x+E*x*x)"
"sin(A+B*x)*exp(C+D*x)+E+F*x"

从这个过程中我希望至少得到常量 (A、B、C 等) 的值，以及关于匹配适合度的统计数据。虽然有商业程序来做这个工作，但是我预期在现代语言库中应该能够找到一些常见的方法来完成期望表达式的拟合。我怀疑 SciPy 的优化部分可能能够做到这一点，但我看不到它让我定义一个方程。同样，在 R 中似乎也找不到我想要的东西。是否存在我正在寻找的内容，还是我需要自己制作？如果它存在并且我只是找不到，那么我很犹豫要自己制作。

编辑: 我想这样做是为了比 LAB Fit 更好地控制整个过程。LAB Fit 的用户界面很差。我还想能够将范围分成多个部分，并使用不同的曲线表示范围的不同部分。最终结果必须能够在速度上击败带有线性插值的查找表，否则我不感兴趣。

在我的当前问题集中，我有三角函数或 exp()，我需要每秒实时执行 352,800 次（并且仅使用 CPU 的一小部分）。所以我绘制曲线并使用数据驱动曲线拟合器，以获得更便宜的近似值。在过去，LUTs 几乎总是解决方案，但现在跳过内存查找并进行近似有时会更快。

- Nosredna

---

你意识到这从统计学角度来看是一个非常糟糕的想法吗？如果你只是想要一个灵活的数据拟合，可以使用像loess、样条或广义加性模型这样的灵活模型。 - hadley

即使将范围分解成更小的范围也是我必须小心处理的成本。我可以访问各种各样的音频数据插值器，但它们通常对我来说计算密集度太高了。一般来说，一旦我开始将范围分解成片段，我最好使用LUT。在DSP应用中，曲线的近似仍然非常有用。 - Nosredna

Answer 1

您的第一个模型实际上在三个参数上是线性的，可以使用R进行拟合。

 fit <- lm(y ~ x + I(x^2), data=X)

这将会为您提供三个参数。

第二个模型也可以使用R中的nls()进行拟合，但需要提供起始值等常规注意事项。在优化中，统计问题并不一定与数值问题相同，无论选择哪种语言，都不能仅仅优化任何函数形式。

- Dirk Eddelbuettel

---

您最好使用y ~ poly(x, 2)或y ~ ns(x, 2)来得到更好的结果。 - hadley

Answer 2

总的来说，回答你关于R中参数估计的问题（不考虑你提到的方程式的具体内容），我认为你需要使用'nls()'或'optim()'函数... 'nls'是我的首选，因为它为每个估计参数提供误差估计，当它失效时我会使用'optim'。如果你有自己的x,y变量：

out <- tryCatch(nls( y ~ A+B*x+C*x*x, data = data.frame(x,y), 
                start = c(A=0,B=1,C=1) ) ,
                error=function(e) 
                optim( c(A=0,B=1,C=1), function(p,x,y)  
                      sum((y-with(as.list(p),A + B*x + C*x^2))^2), x=x, y=y) )

获取系数，可以使用以下方式：

getcoef <- function(x) if(class(x)=="nls") coef(x) else x$par
getcoef(out)

如果您想在“nls”情况下获取标准误差，

summary(out)$parameters

帮助文件和r-help邮件列表中包含了许多关于每个最小化算法的具体实现（默认情况下在上面的每个示例中使用）及其适用于手头方程的具体形式的讨论。某些算法可以处理框约束，而另一个名为constrOptim()的函数将处理一组线性约束。这个网站也可能有所帮助：http://cran.r-project.org/web/views/Optimization.html

- hatmatrix

---

我可以将公式作为字符串输入吗？ - Nosredna

是的 - 用类似 as.formula(paste("y","A+Bx+Cx^2",sep="~")) 的方式即可实现。 - hatmatrix

在nls的情况下，类似于eval(parse(text=sprintf("sum((y-%s)^2)","A+Bx+Cx^2")))这样的优化应该可以工作（sprintf结构被显示出来，以便您可以插入所需的公式）。 - hatmatrix

Answer 3

看看GNU Octave - 在它的polyfit()和非线性约束求解器之间，应该可以构建出适合你问题的东西。

- moonshadow

---

我有时候会用Octave。我会看看能不能弄清楚。 - Nosredna

Answer 4

你可能不会找到一个单一的程序可以灵活地处理你所举例的多项式和有理函数，更不用说它能解析字符串来确定需要拟合的方程类型了。

对于你的第一个例子，最小二乘多项式拟合器是适合的选择。（你可以自己决定使用几次多项式——二次、三次、四次等等）。对于像你的第二个例子这样的有理函数，如果找不到合适的库，你可能需要自己编写代码。此外，请记住，只要你不需要在拟合数据集的极限之外进行外推，就可以使用足够高次数的多项式来近似你的“真实”函数。

正如其他人所指出的那样，还有其他更广义的参数估计算法也可能证明有用。但这些算法并不完全是“即插即用”的：它们通常需要你编写一些辅助程序，并提供模型参数的初始值列表。对于不幸选择的初始参数估计，这些算法可能会发散或陷入局部最小值或最大值。

- Jim Lewis

---

当我使用商业产品时，通常都不知道哪种方法最有效。LAB Fit会尝试几百种方程式，以查看在我指定的范围内哪个最适合拟合数据。 - Nosredna

我没有考虑到这种情况——如果你正在尝试表征数据集的早期阶段，尝试几个函数族（线性、多项式、幂律、周期性等）确实是有意义的，以查看一个好的拟合可能会是什么样子。我会相应地编辑我的答案。 - Jim Lewis

这些算法可能会发散……是的，我猜商业程序在检查所有选择时出现这种情况时会中止运行。当你一次选择一个表达式时，它们确实允许你调整初始值。 - Nosredna

Answer 5

如果您的系数有限制，并且您知道有一种特定类型的函数适合拟合您的数据，而该函数是一个混乱的函数，标准回归方法或其他曲线拟合方法无法使用，那么您是否考虑过遗传算法？

它们不是我的首选，但如果您正在尝试找到第二个函数的系数，则也许GAs会起作用 - 尤其是如果您使用非标准度量来评估最佳拟合。例如，如果您想找到“(A + Bx + Cx ^ 2) /（Dx + Ex ^ 2）”的系数，使得您的函数和数据之间的平方差之和最小并且结果函数的弧长有一些约束条件，那么随机算法可能是一种很好的方法。

一些注意事项：1）随机算法不能保证最佳解决方案，但它们通常会非常接近。2）您必须小心算法的稳定性。

更长的说明是，如果您处于要从某些函数空间中找到最适合您的数据的函数的阶段（例如，您不会将第二个模型强加于您的数据），则遗传编程技术也可能有所帮助。

- B Rivera

---

那是一个有趣的想法。我会考虑一下。显然，这会很慢。商业程序可以在几秒钟内运行数百个方程式。 - Nosredna

是的，随机算法的另一个缺点是速度可能会比较慢。但好处是，可以得到商业程序无法运行的方程形式。通过允许遗传程序搜索函数类（对这些函数进行操作），例如幂函数、指数函数、对数函数、三角函数、概率密度函数/累积分布函数等，可以找到固定方程集合未给出的解决方案。但是，这需要一定的编码工作，可能不值得投入。 - B Rivera

我总是准备好迎接一次唐吉诃德式的冒险。 - Nosredna

Answer 6

在 R 里，这很容易。

内置的方法名叫做 optim()。它需要一个潜在参数的起始向量和一个函数作为参数。你需要构建自己的错误函数，但这非常简单。

然后你可以像这样调用它：out = optim(1, err_fn)

其中 err_fn 是

err_fn = function(A) {
    diff = 0;
    for(i in 1:data_length){
      x = eckses[i];
      y = data[i];
      model_y = A*x;
      diff = diff + ( y - model_y )^2
    }
    return(diff);
}

这里假设您有一个包含x和y值的向量eckses和data。根据需要更改model_y行，甚至添加更多参数。

它可以很好地处理非线性问题，我用它来处理四维e^x曲线，速度非常快。输出数据包括拟合结束时的误差值，这是一个衡量拟合程度的指标，以平方差的总和形式给出（在我的err_fn中）。

编辑： 如果您需要将模型作为字符串输入，则可以让用户界面构建整个模型拟合过程作为R脚本并加载运行。R可以从STDIN或文件中获取文本，因此应该不难构建此函数的字符串等效项，并使其自动运行optim。

- Karl

---

我不使用nls有两个原因，第一，我喜欢手工制作误差函数以优化，第二，我实际上并没有很丰富的R经验。所以nls（）做的恰好是我在上面写的？很棒。 - Karl

我的最终目标是将一个字符串列表传递给代码，并尝试使用它们来找出最佳匹配。 - Nosredna