1

x と y の両方に誤差がある直線の最小二乗フィッティングを実行する以下の関数 (後のブートストラップ用) を高速化しようとしています。主なハングアップは while ループにあると思います。関数の入力値は観測値xおよびであり、これらの値およびy の絶対不確実性です。sxsy

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)
}
4

1 に答える 1

3

いくつかの簡単な変更を行うだけで、関数を高速化できます。主に、そこにある必要のない while ループの外に移動する必要があります。たとえばsolve、同じデータに対して 2 回実行するとします。sfitまた、 while ループの最後の反復でのみ使用する場合は、反復ごとに計算します。

これが私のコードです：

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常にゼロに固定されていることに気付きました。これはもしかして間違いでしょうか？

また、 aを 1 つの列でacbindに変換するためによく使用することにも気付きました。すべてのベクトルは自動的に 1 列の行列と見なされるため、多くの余分なコードを避けることができます。vectormatrix

@joranが述べたように、許容レベルが非常に小さく設定されているため、収束に時間がかかります。公差を大きくすることを検討してください。

于 2012-04-23T18:47:19.777 に答える