使用向量化和矩阵在 R 中回归答案

【问题标题】：Regression in R using vectorization and matrices使用向量化和矩阵在 R 中回归
【发布时间】：2012-01-23 17:16:11
【问题描述】：

我在 R 中有一个使用矩阵的向量化 Q。我有 2 个列需要使用某些索引对每个列进行回归。数据是

matrix_senttoR = [ ...
                  0.11 0.95
                  0.23 0.34
                  0.67 0.54
                  0.65 0.95
                  0.12 0.54
                  0.45 0.43 ] ;
indices_forR = [ ...
            1
            1
            1
            2
            2
            2 ] ;

矩阵中的 Col1 是 MSFT 和 GOOG 的数据（每行 3 行），Col2 是基准 StkIndex 在相应日期的返回值。数据是矩阵格式，因为它是从 Matlab 发送的。

我现在使用

slope <- by(    data.frame(matrix_senttoR),   indices_forR,   FUN=function(x)  
                         {zyp.sen (X1~X2,data=x) $coeff[2] }      ) 
betasFac <- sapply(slope , function(x) x+0)

我正在使用上面的 data.frame，因为我无法使用 cbind()。如果我使用 cbind() 那么 Matlab 会给出一个错误，因为它不理解该数据格式。我正在从 Matlab (http://www.mathworks.com/matlabcentral/fileexchange/5051) 内部运行这些命令。您可以将zyp (zyp.sen) 替换为lm。

BY 这里很慢（可能是因为数据帧？）。有更好的方法吗？ 150k 行数据需要 14 秒以上。我可以在 R 中使用矩阵向量化吗？谢谢。

【问题讨论】：

如果您只是在运行回归，为什么还要将代码从 MATLAB 传递到 R？ Stats 工具箱中的 MATLAB regress 函数可以解决问题。
对代码进行一些分析以查看减速的位置也是一个好主意。您需要知道by 占用了多少时间，建模函数占用了多少时间，以及在 MATLAB 和 R 之间传递数据花费了多少时间。
@Richie -> 这是因为我正在尝试进行非参数回归，特别是使用 zyp 库包。我所有的数据都在 Matlab 中。我唯一的选择是自己在 Matlab 中设计 Theil-Sen Regressor！

标签： r matrix dataframe vectorization regression

【解决方案1】：

这可以很容易地移到评论中，但是：

有几点需要考虑，我倾向于避免使用by() 函数，因为它的返回值是一个时髦的对象。相反，为什么不将您的 indices_forR 向量添加到 data.frame 中？

df <- data.frame(matrix_senttoR) 
df$indices_forR <- indices_forR

plyr 包从这里开始工作：

ddply(df,.(indices_forR),function(x) zyp.sen(X1~X2,data=x)$coeff[2])

您可以使用 doMC 或 doSnow 以及 .parallel=TRUE 到 ddply 的参数轻松地 multi-thread 此操作。

如果速度是目标，我也会学习data.table 包（它包装了 data.frame 并且速度更快）。另外，我假设慢的部分是zyp.sen() 调用而不是by() 调用。在多个内核上执行将加快这一进程。

> dput(df)
structure(list(X1 = c(0.11, 0.23, 0.67, 0.65, 0.12, 0.45), X2 = c(0.95, 
0.34, 0.54, 0.95, 0.54, 0.43), indices_forR = c(1, 1, 1, 2, 2, 
2)), .Names = c("X1", "X2", "indices_forR"), row.names = c(NA, 
-6L), class = "data.frame")

> ddply(df,.(indices),function(x) lm(X1~X2,data=x)$coeff[2])
  indices         X2
1       1 -0.3702172
2       2  0.6324900

【讨论】：

-> 步骤 ddply 在像 evalR('slope
@Maddy 听起来像是 Matlab 错误而不是 R 错误。不确定 zyp.sen() 来自哪个包，但在 R 本身中使用 lm() 效果很好。
-> 谢谢贾斯汀。但我想我现在坚持使用 zyp。它是非参数的，我必须专门使用它。我知道这是基于 matlab 的错误。我不是很精通 R，所以把这个 Q 希望有一个解决方案。

【解决方案2】：

我仍然认为您从 MATLAB 迁移到 R 再返回会使事情变得过于复杂。传递 150k 行数据肯定会大大减慢速度。

zyp.sen 移植到 MATLAB 实际上非常简单。给你：

function [intercept, slope, intercepts, slopes, rank, residuals] = ZypSen(x, y)
% Computes a Thiel-Sen estimate of slope for a vector of data.

n = length(x);

slopes = arrayfun(@(i) ZypSlopediff(i, x, y, n), 1:(n - 1), ...
    'UniformOutput', false);
slopes =  [slopes{:}];
sni = isfinite(slopes);
slope = median(slopes(sni));

intercepts = y - slope * x;
intercept = median(intercepts);

rank = 2;
residuals = x - slope * y + intercept;

end


function z = ZypSlopediff(i, x, y, n)

z = (y(1:(n - i)) - y((i + 1):n)) ./ ...
    (x(1:(n - i)) - x((i + 1):n));

end

我使用 R 的example(zyp.sen) 进行了检查，它给出了相同的答案。

x = [0 1 2 4 5]
y = [6 4 1 8 7]
[int, sl, ints, sls, ra, res] = ZypSen(x, y)

您确实应该做一些进一步的检查，以确保。

【讨论】：