0

我试图找出一种更有效的方法来分析具有 4 个自变量（x1、x2、x3 和 x4）和 1 个响应变量（y，可以是 0 到 1 之间的任何数字）的数据集。我正在尝试将我的数据拟合到具有 6 个参数（3 个斜率 [m1、m2 和 m3] 和 3 个截距 [b1、b2 和 b3]）的模型中。

每个点的 y 响应通过求解 y 来计算，其中以下等式等于 0：

-x4+(x1/((log(-y/(y-1))-b1)/m1))
          +(x2/((log(-y/(y-1))-b2)/m2))
          +(x3/((log(-y/(y-1))-b3)/m3))

据我所知（如果有更好的方法，请纠正我），实现这一点的最佳方法是用 最小化上述方程的绝对值optimize()。这是一个具有任意参数和 x 值的可重现示例：

#model parameters
b1=-8
b2=-10
b3=-15
m1=2
m2=25
m3=50

#independent variables
x1=3.9
x2=.02
x3=.01
x4=1

a=function(y) abs(-x4+(x1/((log(-y/(y-1))-b1)/m1))
      +(x2/((log(-y/(y-1))-b2)/m2))
      +(x3/((log(-y/(y-1))-b3)/m3)))

y=optimize(a,c(0,1))

使用这些输入，y 评估为 ~0.617。

很简单，但我必须对 500 个数据点中的每一个都执行此操作（每个数据点都有唯一的 x1、x2、x3、x4 组合，但都具有相同的参数 b1、b2、b3、m1、m2 和 m3）。我目前正在使用vapply()，但似乎必须有一种更有效的方法。但是，我不知道向量化优化问题的方法。

如果它在这里结束，它不会那么糟糕（所有 500 分都在不到一秒的时间内评估vapply()）optimization()。但这仅计算一组给定参数的 y 变量。当我尝试使用DEoptim()最大化对数似然来优化参数 b1、b2、b3、m1、m2 和 m3 时出现问题（然后我使用参数 fromDEoptim()作为参数微调的起始参数optim()）。不用说，由于每次参数迭代都需要评估 500 个优化问题，这需要一段时间。奇怪的是，当我手动运行 vapply 以优化一组参数（即，相当于 DEoptim 的单次迭代）时，它花费的时间不到一秒，但每次迭代DEoptim()（这也应该花费不到一秒钟）大约需要 10 秒钟。我不知道为什么vapply()要花这么长时间（有什么想法吗？），但我希望有一种方法可以更有效地解决所有 500 个方程。

在诉诸vapply()and之前，我尝试了一些替代方法，包括使用(而不是)optimize()迭代计算 y （同时在所有 500 个点上），它有点快（每次迭代约 6 秒），但我认为我应该每次迭代能够实现 <1 秒（就像我手动运行一个实例时一样）。如果找不到更好的选择，我可能会回到.Reduce()optimize()vapply()Reduce()

非常感谢任何帮助或见解。谢谢！

编辑：为了澄清，指定的方程在每个点为指定的参数和 x 值（总共 500 个数据点，由 x1、x2、x3 和 x4 的唯一组合组成，所以有 500 个唯一的 y 值要计算；y 是每个实例中唯一的未知数）。总体目标是通过最大似然优化参数以获得最适合观察到的 y 值的预测 y 值。在参数优化的每次迭代中，都必须重新计算 500 y 值，因为 m1 m2 m3 和 b1 b2 b3 参数已经改变。

4

1 回答 1

3

最好使用平方而不是 abs() 作为目标函数，因为这样一阶导数是连续的。你可以通过这样做看到这一点

curve(a,.001,.999)

并通过省略 abs 并定义此功能

b=function(y) -x4+(x1/((log(-y/(y-1))-b1)/m1))+
              (x2/((log(-y/(y-1))-b2)/m2))+
              (x3/((log(-y/(y-1))-b3)/m3))

并绘制这个函数

curve(b,0.001,.9999)

一般来说，使用优化算法找到方程组的解并不是一个好主意，因为它会寻找任何最小值（全局和局部）。你想要一个最小值为 0。

所以最好使用非线性方程求解器。有一个包nleqslv可以做到这一点（注意：我是那个包的作者）。

Vectorize由于您的函数已经矢量化，因此无需使用或其他矢量化方法。

定义一个函数（与你的没有f相同，但效率更高）aabs

f <- function(y) {
    tmp <- log(-y/(y-1))
    -x4+(x1/((tmp-b1)/m1))+(x2/((tmp-b2)/m2))+(x3/((tmp-b3)/m3))

}

并定义一个计算雅可比的函数

fjac <- function(x) { h <- 0.00001*x; diag((f(x+h)-f(x))/h) }

这可以非常有效地完成，因为返回值的每个元素f仅取决于输入向量的相应元素y。

y对于每个参数配置，数据向量和解决方案的起始值y可以通过下式计算

z <- nleqslv(ystart,f,fjac, method="Newton")
y <- z$x

您必须使用方法Newton，因为方法Broyden在这种情况下不起作用。你可以试试这个例子

K <- 500
x1 <- x1 + c(0,runif(K-1,.1*x1,.3*x1))
x2 <- x2 + c(0,runif(K-1,.01*x2,.03*x2))
x3 <- x3 + c(0,runif(K-1,.01*x3,.03*x3))
x4 <- x4 + c(0,runif(K-1,.01*x4,.03*x4))
nleqslv(rep(.3,K),f,fjac, method="Newton")

在我的电脑上，这大约需要 0.08 秒。

于 2013-08-25T09:36:57.643 回答