多元多项式的拟合分数答案

【问题标题】：Fitting fractions of multivariate polynomial多元多项式的拟合分数
【发布时间】：2017-02-03 17:31:03
【问题描述】：

如何拟合以下形式的多元分数多项式：

给定一个函数：y = f(x,z)，两个变量x 和z 的函数。更具体地说，它的形式是：

y = (x^2 + x^3)/(z^2 + z^3)

分子是预测变量x的3次多项式，分母也是某个预测变量z的3次多项式。

我想为每个预测变量 x 和 z 拟合多项式，即我需要找到系数 A、B、C、D：

y = (A*x^2 + B*x^3)/(C*z^2 + D*z^3)

基本上，y 是两个 3 次多项式的比值。如何拟合这样的函数？

下面的数据框示例。我无法发布完整的数据框，因为它有超过 1000 行。

y = c(-4.10594369806545, -4.23691458506868, -4.24690667936422, -3.53677470254628, -4.30406509320417, -4.19442802077908, -4.66857169733859, -2.82271310942235, -4.19720194766181, 3.52164353473802, -4.3917019001973, -5.41654474791269, 2.87471821731616, -3.85922481986118, -4.25370811223789, -3.57887855889961, -5.33913936106829, -4.11775265312012, -2.89958841300109, -4.18661983833127)

x = c(8.06526520889773, 9.39897529082673,9.07348918922699,7.5522372875608, 9.17294998275762,5.77455154554441, 9.2005930205213, 8.07309119969315, 7.42177579364465,8.18896686364888, 8.07868822922987, 8.50956416425175,9.71269017726113, 7.98378106897745, 7.69893619981345, 8.49576524400262, 8.02224091680654,8.25400859056484, 7.58171964012531, 8.35655484545343)

z = c(2.56494935746154, 4.99043258677874, 4.43081679884331,3.66356164612965,4.53259949315326,1.79175946922805,4.23410650459726, 5.38449506278909,3.13549421592915,4.34380542185368, 3.43398720448515,2.77258872223978,6.94985645500077,3.97029191355212, 3.40119738166216,4.39444915467244,2.19722457733622,3.91202300542815,4.06044301054642, 3.87120101090789)

dat = data.frame(cbind(y=y,x=x,z=z))

更新：

致电nls:

nls(y~(a*(x**2) + b*(x**3))/(c*(z**2) + d*(z**3)), dat, start=list(a=1,b=1,c=1,d=1))

【问题讨论】：

使用超过 10^3 行的扩展数据集，我看到了我刚刚发布的错误。使用我在此处提供的示例，我也遇到了渐变问题。
感谢您的尝试。我仍然得到 y 和 z 的参数初始化错误。
你是否经历过参数的迭代？你是怎么得到 C=10 的？
但是，我可以通过浏览参数空间并计算 R^2 和其他模型优度指标来找到最佳模型。你认为这样可以找到正确的参数吗？
我认为另外使用optim() 会更好。 abc.f <- function(par) sum(((par[1]*x**2 + par[2]*x**3) / (par[3]*z**2 + z**3) - y)^2); optim(par=c(1, 1, -7), fn = abc.f, control=list(maxit = 1e+4, reltol=1e-12))（示例）

标签： r regression nls non-linear-regression

【解决方案1】：

这是一个很好的问题。至少我从中学到了一些东西。但是，我有一种感觉，这个问题更多地围绕这个特定任务（大学作业？）的解决方案，而不是一个一般性问题。

但让我分享一种解决方法：我们这里有什么

可以简化为

求解 y^\theta 在数值上变得更易于管理。正如我们所看到的（在非常努力地尝试并未能解决非线性问题之后）它实际上是两个线性模型的划分。所以一种方法是分别估计这两个问题的系数。也就是说，我们固定系数a 和b 找到c 和d，然后使用c 和d 找到a 和b。

以下代码求解系数

首先是一些数据

library(dplyr)

sampleData <- data.frame(x = runif(100, -100, 100), z = runif(100, -100, 100)) %>%
  mutate(y = ( (-2 * x^2) + (5 * x^3) ) /  (-4 * z^2 + 6 * z^3)) %>%
  mutate(zxfactor = z^2/x^2,
         yy = y * zxfactor)

现在我们求解yy。有一些随机的起始值...

init2 <- structure(runif(4, -10, 10), names=c("A", "B", "C", "D"))
coefab <- init2[c("A", "B")]
coefcd <- init2[c("C", "D")]

...我们需要为a 和b 拟合线性模型

以及c 和d 的线性模型

# don't use for loop but determine a terminal condition... but i'm too lazy :-)
for(i in 1:100) {
  # make linear prediction using coeff. c and d
  sampleData <- sampleData %>%
    mutate(yab = yy * (coefcd[1] + coefcd[2] * z))
  # and fit a model for a and b
  coefab <- coef(lm(yab ~ x, sampleData))
  # then make a linear prediction using coeff. a and b
  sampleData <- sampleData %>%
    mutate(ycd = (coefab[1] + coefab[2] * x) / yy)
  # and fit a model for c and d
  coefcd <- coef(lm(ycd ~ z, sampleData))
} # repeat until satisfied

coefab
coefcd

我们对找到的系数满意吗？让我们检查一下：

optimFun <- function(params, out, x, z) {
  res <- (params[1] + params[2]*x)/(params[3] + params[4]*z)
  return( sqrt(sum( (out - res)^2 )) )
}

optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy)

> optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy)
[1] 1.951043e-12

事实上，我们是因为建模函数估计和数据yy（按比例缩放的）之间的差异接近于零。由于问题是超定的，因此不同的迭代会导致不同的参数估计。（也许有人可以更详细地解释一下）

评论：

估计值比nls 确定的估计值更接近于零
如果您想使用optim，您可以使用optimFun。为了更快地收敛，您甚至可以定义导数函数
这是一个很好的问题，它显示了一般优化可能会失败的地方，并且考虑手头的问题总是值得的。
尝试lm(yyz ~ x, data = sampleData %>% mutate(yyz = yy*(-4 + 6*z)))，它将返回a = -2 和b = 5 的精确值。给定任意两个参数，您可以找到最小化函数的匹配对。

【讨论】：