目的

観察値と予測値の差の絶対値の p 乗の和を最小にするという方法による回帰直線のパラメータを推定する。
p が 2 のときは，普通の最小二乗法である。
p が 1 のときは，絶対値の和を最小とするものである。

使用法

least.sum.abs(x, y, n=1, p=1, control=list())

# print メソッド
print.least.sum.abs(obj, digits=5, sig=0.95)

# plot メソッド
plot.least.sum.abs(obj, posx="topleft", posy=NULL, xlab=obj$names.xy[1], ylab=obj$names.xy[2],	hist=FALSE, ...)

引数

x             独立変数ベクトル
y             従属変数ベクトル
n             ブートストラップ法で信頼区間を求めるときの回数
              デフォルトは 1。この場合には信頼区間を求めない
p             デフォルトは 1。2 のときは，普通の最小二乗法
obj           "least.sum.abs" オブジェクト
digits        表示桁数
sig           ブートストラップ法による信頼区間を求めるときの信頼度（デフォルトで 95% 信頼限界）
posx, posy    legend 関数のための位置引数（"topleft", "bottomright" なども可）
xlab, ylab    軸の名前
hist          ブートストラップ法による結果をヒストグラム表示するかどうか
...           plot 関数に渡されるその他の引数

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/least.sum.abs.R", encoding="euc-jp")

# 観察値と予測値の差の絶対値の p 乗の和を最小にするという方法による回帰直線のパラメータ推定
least.sum.abs <- function(   x,                                      # 独立変数ベクトル
                                y,                                      # 従属変数ベクトル
                                n=1,                                    # ブートストラップ法で信頼区間を求めるときの回数
                                p=1,                                    # 1 のとき絶対値の和を最小にする。2 のときは普通の最小二乗法
                                control=list())                         # optim 関数への引数
{
        least.sum.abs0 <- function(x, y)                             # 1 組のデータについて，切片と傾きの推定値を計算する
        {
                evaluate.error <- function(par)
                {
                        return(sum(abs(y-par[2]*x-par[1])^p))
                }
                if (p < 1) {
                        warning("p は 1 以上の値でなければなりません。p=1 として実行します。")
                        p <- 1
                }
                ans <- optim(c(0, 0), evaluate.error, control=control)
                return(c(Intercept=ans$par[1], Slope=ans$par[2]))
        }
        Driver <- function(x, y)                                     # ブートストラップ法のためのドライバー
        {
                n <- length(x)
                suffix <- sample(n, n, replace=TRUE)                 # リサンプリング
                return(least.sum.abs0(x[suffix], y[suffix]))            # リサンプリングしたデータについてパラメータを推定
        }
        names.xy <- c(deparse(substitute(x)), deparse(substitute(y)))        # 変数名を控えておく
        OK <- complete.cases(x, y)                                   # 欠損値を持つケースを除く
        x <- x[OK]
        y <- y[OK]
        ans <- list(coefficients=least.sum.abs0(x, y),                       # 引数に対してパラメータを推定する
                    names.xy=names.xy, x=x, y=y)
        if (n > 1) {
                ans2 <- replicate(n, Driver(x, y))                   # ブートストラップを n 回実行
                ans <- append(ans, list(intercepts=ans2[1,], slopes=ans2[2,]))
        }
        class(ans) <- "least.sum.abs"                                        # print, plot メソッドのためにクラス名をつけておく
        return(ans)
}
# print メソッド
print.least.sum.abs <- function(obj,                                 # "least.sum.abs" オブジェクト
                                digits=5,                               # 表示桁数
                                sig=0.95)                               # 信頼度
{
        if (length(obj) == 4) {
                cat("Intercept:", round(obj$coefficients[1], digits),
                    "    Slope:", round(obj$coefficients[2], digits), "\n")
        }
        else {
                alpha <- (1-sig)/2
                LCL <- c(quantile(obj$intercepts, alpha), quantile(obj$slopes, alpha))
                UCL <- c(quantile(obj$intercepts, 1-alpha), quantile(obj$slopes, 1-alpha))
                ans <- data.frame(obj$coefficients, LCL, UCL)
                dimnames(ans) <- list(c("Intercept:", "Slope:"),
                                      c("Estimate", paste(c(alpha, 1-alpha), "%", sep="")))
                print(round(ans, digits=digits))
        }
}
# plot メソッド
plot.least.sum.abs <- function(obj,                                  # "least.sum.abs" オブジェクト
                        posx="topleft", posy=NULL,                      # legend 関数のための位置引数
                        xlab=obj$names.xy[1], ylab=obj$names.xy[2],     # 軸の名前
                        hist=FALSE,                                     # ヒストグラムを描くとき TRUE にする
                        ...)                                            # その他の任意の plot 関数の引数
{
        if (hist && length(obj) == 6) {                                   # ブートストラップの結果を，hist=TRUE のときに，ヒストグラムで表示する
                layout(matrix(1:2, 2))
                hist(obj$intercepts, xlab="Intercept", main="", right=FALSE)
                hist(obj$slopes, xlab="Slope", main="", right=FALSE)
                layout(1)
        }
        else {                                                          # 散布図と least.sum.abs 法の回帰直線と直線回帰式を表示する
                plot(obj$x, obj$y, xlab=xlab, ylab=ylab, ...)
                abline(obj$coefficients)
                abline(lm(obj$y~obj$x), lty=2, col=2)
                legend(posx, posy, legend=c("least sum abs", "linear regression"), lty=1:2, col=1:2)
        }
}


使用例

> set.seed(123)

> d <- gendat2(100, 0.7)*15+50

> (a <- least.sum.abs(d[,1], d[,2])) # 点推定値のみを求める

Intercept: 16.07177     Slope: 0.66218 

> plot(a, xlab="x", ylab="y", pch=19, xlim=c(0,100),ylim=c(0,100)) # 散布図と回帰直線



> (b <- least.sum.abs(d[,1], d[,2], n=1000)) # 点推定値と信頼区間を求める

           Estimate  0.025%   0.975%
Intercept: 16.07177 6.14666 29.87396
Slope:      0.66218 0.41653  0.91720

> plot(b, hist=TRUE) # ブートストラップパラメータの分布図