【问题标题】:QR decomposition different in lm and biglm?lm和biglm中的QR分解不同?
【发布时间】:2012-10-26 06:48:24
【问题描述】:

我正在尝试从 biglm 中使用的 QR 分解中恢复 R 矩阵。为此,我使用了 vcov.biglm 中的一部分代码并将其放入如下函数中:

qr.R.biglm <- function (object, ...) {
    # Return the qr.R matrix from a biglm object
    object$qr <- .Call("singcheckQR", object$qr)
    p <- length(object$qr$D)
    R <- diag(p)
    R[row(R) > col(R)] <- object$qr$rbar
    R <- t(R)
    R <- sqrt(object$qr$D) * R
    dimnames(R) <- list(object$names, object$names)
    return(R)
}

更具体地说,我试图从基本包中获得与使用 qr.R 相同的结果,该包用于类“qr”的 QR 分解,例如包含在 lm 类 (lm$qr) 中的那些。基函数代码如下:

qr.R <- function (qr, complete = FALSE) {
    if (!is.qr(qr)) 
        stop("argument is not a QR decomposition")
    R <- qr$qr
    if (!complete) 
        R <- R[seq.int(min(dim(R))), , drop = FALSE]
    R[row(R) > col(R)] <- 0
    R
}

我设法从样本回归中得到相同的结果,除了符号。

x <- as.data.frame(matrix(rnorm(100 * 10), 100, 10))
y <- seq.int(1, 100)
fit.lm <- lm("y ~ .", data =  cbind(y, x))
R.lm <- qr.R(fit.lm$qr)

library(biglm)
fmla <- as.formula(paste("y ~ ", paste(colnames(x), collapse = "+")))
fit.biglm <- biglm(fmla, data = cbind(y, x))
R.biglm <- qr.R.biglm(fit.biglm)

比较两者,很明显绝对值匹配,但符号不匹配。

mean(abs(R.lm) - abs(R.biglm) < 1e-6)
[1] 1
mean(R.lm - R.biglm < 1e-6)
[1] 0.9338843

我不太明白这是为什么。我希望能够从 biglm 获得与 lm 相同的 R 矩阵结果。

【问题讨论】:

    标签: r matrix lm qr-decomposition


    【解决方案1】:

    两个 R 矩阵之间的区别在于,biglm 显然执行它的旋转使得 R 的对角线元素都是正数,而 lm(或者,实际上,它的例程call) 没有施加这样的约束。 (一种策略或另一种策略应该没有数值优势,所以区别只是惯例之一,AFAIKT。)

    您可以通过自己施加附加约束来使 lm 的结果与 biglm 的结果相同。我会使用将列乘以 1 或 -1 的反射矩阵,这样对角线元素最终都是正数:

    ## Apply the necessary reflections
    R.lm2 <- diag(sign(diag(R.lm))) %*% R.lm
    
    ## Show that they did the job
    mean(R.lm2 - R.biglm < 1e-6)
    # [1] 1
    

    【讨论】:

    • 谢谢乔希。这很有用。但是,我正在尝试解决相反的问题:如何使 biglm R 矩阵与 lm R 矩阵匹配?如果可能的话。
    • @user3327 -- 好的。我想我误解了你帖子的最后一句话。我将无法再提供任何帮助,但如果其他人有兴趣,我会留下来......只是出于好奇,你为什么关心两个 R 矩阵的列的符号是否匹配?
    • 谢谢乔希。只是为了提供更多背景知识:我正在尝试在一些现有代码中将lm 替换为biglm,因为我有一个大型数据集,我想对其执行回归并且将无法使用lm .我的程序进一步使用了这个 R 矩阵,如果我可以简单地提取相同的矩阵,而不是必须想出一个新的公式,那就太好了。如果确实如此,那就这样吧:-)我只是想先问问。
    猜你喜欢
    • 2014-03-11
    • 2016-12-09
    • 2013-11-01
    • 2014-06-09
    • 1970-01-01
    • 2021-06-06
    • 2012-06-13
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多