R：加速多个 lm()答案

【问题标题】：R: Speed up multiple lm()R：加速多个 lm()
【发布时间】：2016-04-27 01:55:10
【问题描述】：

我想估计一个非线性模型的参数。

模型方程为Z = A * exp(- a * X) + B * exp(- b * Y) + C

X 和 Y 是预测变量
A、B、a、b 是要估计的参数

我所做的是在进行线性回归之前通过指数转换将模型转换为线性问题：

对于介于 0 和 1 之间的 a 和 b，我计算 exp_x = exp(- a * X) 和 exp_y = exp(- b * Y)
我做了一个线性回归Z ~ exp_x + exp_y

正如我们在这个模拟中看到的那样，它工作得非常好

x = 1:10
y = 1:10

combination = expand.grid(x = x, y = y)

df = data.frame(
  X = combination$x,
  Y = combination$y,
  Z = 2 * exp(-0.3 * combination$x) + 
      5 * exp(-0.6 * combination$y) + 
      rnorm(n = 100, mean = 0, sd = 0.1 )
)

a_hat = 0
b_hat = 0
best_ols = NULL
best_rsquared = 0

for (a in seq(0.01, 1, 0.01)){
  for (b in seq(0.01, 1, 0.01)){  

    df$exp_x = exp(- a * df$X)
    df$exp_y = exp(- b *df$Y)

    ols = lm(data = df, formula =  Z ~ exp_x + exp_y)
    r_squared = summary(ols)$r.squared

    if (r_squared > best_rsquared){
      best_rsquared = r_squared 
      a_hat = a
      b_hat = b
      best_ols = ols
    }    
  }
}

a_hat 
b_hat 
best_ols
best_rsquared 

> a_hat 
[1] 0.34
> b_hat 
[1] 0.63
> best_ols

Call:
lm(formula = Z ~ exp_x + exp_y, data = df)

Coefficients:
(Intercept)        exp_x        exp_y  
     0.0686       2.0550       5.1189  

> best_rsquared
[1] 0.9898669

问题：这很慢

大约需要 10 秒，我需要在其他数据帧上执行数千次。

我怎样才能大大加快速度？

【问题讨论】：

标签： r performance loops lm

【解决方案1】：

也许改用nls。由于您没有set.seed()，因此无法看到我们的预测是否相似，但至少我在编辑后得到了a 和b 的“正确”估计：

nmod <- nls( Z ~ A*exp(-a*X)+B*exp(-b*Y), data=df, start=list(A=0.5, B=0.5, a=.1,b=.1))

> coef(nmod)
        A         B         a         b 
2.0005670 4.9541553 0.2951589 0.5937909 
#--------
> nmod
Nonlinear regression model
  model: Z ~ A * exp(-a * X) + B * exp(-b * Y)
   data: df
     A      B      a      b 
2.0006 4.9542 0.2952 0.5938 
 residual sum-of-squares: 0.9114

Number of iterations to convergence: 9 
Achieved convergence tolerance: 5.394e-06

比您的 10 秒体验快得多。这是在一台使用了 8 年的机器上。

> system.time( nmod <- nls( Z ~ A*exp(-a*X)+B*exp(-b*Y), data=df, start=list(A=0.5, B=0.5, a=.1,b=.1)) )
   user  system elapsed 
  0.036   0.002   0.033

【讨论】：