R更快的方法while循環最小二乘函數

Question

我試圖加快下面的函數（用於以後的自舉），它執行最小二乘法擬合一個在x和y中都有錯誤的直線。我認爲主要的掛斷是在while循環中。該函數的輸入值是觀察值x和y以及這些值的絕對不確定性值sx和sy。

york <- function(x, y, sx, sy){ 

    x <- cbind(x) 
    y <- cbind(y) 

    # initial least squares regression estimation 
    fit <- lm(y ~ x) 
    a1 <- as.numeric(fit$coefficients[1]) # intercept 
    b1 <- as.numeric(fit$coefficients[2]) # slope 
    e1 <- cbind(as.numeric(fit$residuals)) # residuals 
    theta.fit <- rbind(a1, b1) 

    # constants 
    rho.xy <- 0  # correlation between x and y 

    # initialize york regression 
    X <- cbind(1, x) 
    a <- a1 
    b <- b1 
    tol <- 1e-15 # tolerance 
    d <- tol 
    i = 0 

    # york regression 
    while (d > tol || d == tol){ 
     i <- i + 1 
     a2 <- a 
     b2 <- b 
     theta2 <- rbind(a2, b2) 
     e <- y - X %*% theta2 
     w <- 1/sqrt((sy^2) + (b2^2 * sx^2) - (2 * b2 * sx * sy * rho.xy)) 
     W <- diag(w) 
     theta <- solve(t(X) %*% (W %*% W) %*% X) %*% t(X) %*% (W %*% W) %*% y 

     a <- theta[1] 
     b <- theta[2] 

     mswd <- (t(e) %*% (W%*%W) %*% e)/(length(x) - 2) 
     sfit <- sqrt(mswd) 
     Vo <- solve(t(X) %*% (W %*% W) %*% X) 
     dif <- b - b2 
     d <- abs(dif) 
     } 

    # format results to data.frame 
    th <- data.frame(a, b) 
    names(th) <- c("intercept", "slope") 
    ft <- data.frame(mswd, sfit) 
    names(ft) <- c("mswd", "sfit") 
    df <- data.frame(x, y, sx, sy, as.vector(e), diag(W)) 
    names(df) <- c("x", "y", "sx", "sy", "e", "W") 

    # store output results 
    list(coefficients = th, 
     vcov = Vo, 
     fit = ft, 
     df = df) 
} 

Answer 1

您的功能可以通過一些簡單的更改加快。首先，你應該將任何東西移出while循環，而不需要在那裏。例如，您在相同的數據上運行兩次solve。此外，在每次迭代中，只有在while循環的最後一次迭代中使用它時，才計算sfit。

這裏是我的代碼：

york.fast <- function(x, y, sx, sy, tol=1e-15){ 
    # initial least squares regression estimation 
    fit <- lm(y ~ x) 
    theta <- fit$coefficients 
    # initialize york regression 
    X <- cbind(1, x) 
    d <- tol 
    # york regression 
    while (d >= tol){ 
     b2 <- theta[2] 
     # w <- 1/sqrt((sy^2) + (b2^2 * sx^2) - (2 * b2 * sx * sy * rho.xy)) # rho.xy is always zero! 
     w <- 1/sqrt(sy^2 + (b2^2 * sx^2)) # rho.xy is always zero! 
     # W <- diag(w) 
     # w2 <- W %*% W 
     w2 <- diag(w^2) # As suggested in the comments. 
     base <- crossprod(X,w2) 
     Vo <- solve(base %*% X) 
     theta <- Vo %*% base %*% y 
     d <- abs(theta[2] - b2) 
    } 
    e <- y - X %*% theta 
    mswd <- (crossprod(e,w2) %*% e)/(length(x) - 2) 
    sfit <- sqrt(mswd) 

    # format results to data.frame 
    th <- data.frame(intercept=theta[1], slope=theta[2]) 
    ft <- data.frame(mswd=mswd, sfit=sfit) 
    df <- data.frame(x=x, y=y, sx=sx, sy=sy, e=as.vector(e), W=diag(diag(w))) 

    # store output results 
    list(coefficients = th, vcov = Vo, fit = ft, df = df) 
} 

一個小測試：

n=225 
set.seed(1) 
x=rnorm(n) 
y=rnorm(n) 
sx=rnorm(n) 
sy=rnorm(n) 

system.time(test<-york.fast(x,y,sx,sy)) # 0.37 s 
system.time(gold<-york(x,y,sx,sy)) # 1.28 s 

我注意到，rho.xy始終固定在零。這可能是一個錯誤？

我注意到你也經常使用cbind將vector轉換成matrix與一列。所有向量都被自動視爲具有一列的矩陣，所以可以避免大量額外的代碼。

正如@joran所提到的，公差水平設置得很小以至於需要很長時間來收斂;考慮使用更大的容差。