＜論 文＞
最小絶対偏差推定量のための効率的算法
近 藤 康 之
1. は じ め に
本稿の目的は，線形回帰モデルにおける最小絶 対偏差推定量（least absolute deviations estimator, LADE）（Amemiya［2, 4.6］）を 計 算 す る た めの算法の効率性について，ブートストラップ法 を用いる場 合 も 含 め て 検 討 す る こ と で あ る。
LADE を求めるための算法としては，線形計画 問題（linear program,LP）によるものが標準的 であり，そこでは単体法（simplex method）に 基づく算法が使用されることが多い（LP と単体 法などの算法については，例えば今野［1］を参 照していただきたい）。これまでに提案された典 型的な方法は，LADE を定義する最小化問題を LP に書き換えた上で，一般的な LP のための単 体法を LADE の特殊な構造を利用して改善した ものである。例えば Armstrong and Kung［3］
は，１回のピボット演算で複数の基底変数と非基 底変数を入れ替えたり，ピボット演算で非基底変 数から基底変数に変わる候補を逐次的に一部の変 数に限定したりすることで，主単体法に基づく有 界変数法を改善している。
LADE の漸近分布は撹乱項の密度関数に依存 するため，推定量の分散を推定することが容易で はない。そのため，ブートストラップ法を用いる ことが効果的である（Greene［5, 16.3.2］）。
単に LADE の値を求めるだけならば LP を１回 解けば済むが，分散も併せて推定するためには，
LP を繰り返し解く必要がある。したがって，よ り効率的な算法の開発が求められる。観察データ
に基づく１回の LP と，ブートストラップ標本に 基づく LP とは構造がよく似ているため，この類 似点を利用することで，効率的な算法を開発し得 ることが予想される。このような観点から，本稿 では，LADE のための特別な算法を開発するよ りも，広く普及している算法（主単体法，双対単 体法，有界変数法）を比較することで，それらの 特徴を明らかにすることを試みる。
数値計算（シミュレーション）に基づく比較に より，主単体法に基づく有界変数法よりも双対単 体法に基づく有界変数法の方が効率的であること，
主単体法に基づく算法の非効率さは初期可能解の 与え方に大きく依存すること，が明らかになった。
2. 最小絶対偏差推定量を定義する線形計画 問題
次の回帰モデルを考える。
＝ β＋ε ＝1, , .
ここで， はスカラの被説明変数， は ×1の 説明変数ベクトル，βは ×1のパラメタベクト ル，εはスカラの撹乱項である。また，上付き 添字 は行列・ベクトルの転置をあらわす。パラ メタ βの最小絶対偏差推定量（LADE） は，
minimize ∑ − with respect to ,
⑴ の解として定義される。ここで は∑ ＝ を満たすウェイトであり，LADE では ＝ ＝
…＝ ＝1である。
Charnes, Cooper and Ferguson［4］が示した
＊ 早稲田大学政治経済学術院
ように，上の最小化問題 ⑴の最適解 は，次の 線形計画問題（LP）の最適基底解 , , の 一部として求められる。
minimize ∑ ＋∑
subject to − ＝ − , 0, 0
＝1, , with respect to , , , , , , .
⑵
LP ⑵を行列・ベクトルを用いて書き直すと，
minimize subject to
＋
＋ − ＝ , 0 , 0 ,
⑶ となる。また，この双対問題は，
maximize subject to
ξ ξ＝0 ,
− ξ ,
である。さらに，この問題を標準的な有界変数問 題（すべての変数の下限はゼロで共通）のための 算法を適用しやすいように，ζ＝ξ＋ として，
maximize subject to
ζ
ζ＝ , 0 ζ 2 ,
⑷ と書き換えておく（目的関数における定数部分
− を省略した）。
主問題⑶が ＋2 変数 制約式であるのに対 して，双対問題⑷は 変数 制約式である。単体 法に基づく算法の効率性はピボット演算の効率性 に大きく依存し，制約式が多いほど１回のピボッ ト演算の計算量は多くなる。したがって，上の２ つの LP を（非負変数ではなく有界変数であるこ とを適切に考慮した）単体法によって解く限りに おいては，制約式の少ない双対問題⑷の方が計算 量が少ないと言える。以下では，双対問題⑷のみ を考察の対象とする。
3. 単体法に基づく算法と初期可能基底
有界変数を持つ LP のための単体法として，主 単体法に基づく算法（今野［1, 9.1］）と双対 単体法に基づく算法（Maros［6, pp.5−6］）を 比較の対象とする。主単体法に基づく算法のため
には，初期基底として主可能基底を求める必要が あるが，これは容易には求められない。したがっ て，２段階単体法と同様の方法で初期可能基底を 求める（今野［1, pp.141‑142］）。他 方，双 対 単 体法に基づく算法のための初期基底としては，次 に述べるように，比較的容易に双対可能基底を求 めることができるので，これを用いる。
のある基底集合を ⊂ ＝ 1, , とする。
す な わ ち， ＝ で あ り， 個 の 列 ベ ク ト ル
∈ を並べた基底行列 は非特異である。
以下では，文脈から判断して誤解の生じない限り，
基底集合 と基底行列 の両方を基底と呼ぶ。
非基底 ＝ の − 個の列を並べた ×
− 行列を とする。 の基底 と非基底 への分割に対応して，
＝ , ＝ ,
ζ＝ ζζ , ＝ と分割する。
LP ⑷の制約式 ζ＝ と目的関数 ζは，
この分割に沿って次のように書き換えられる。
ζ＋ ζ＝ ＋
ζ＝ ＋ − ζ
ζ＝ ＋ −ζ ,
ζ＝ ζ＋ ζ
＝ ＋ −ζ ＋ ζ
＝ ＋
＋ − ζ.
また，被約費用ベクトル を，
＝ − ⑸
と定義する。基底 に対する基底解は，
ζ＝ ＋ −ζ ,
ζ＝要素毎に０または2 , ⑹ とあらわされる。非基底変数ζについて，
ζ＝ 0 2
0のとき
＞0のとき ∈ とすると，この基底解は双対可能である。
対応する LP ⑶の基底解 , , は，残差ベ クトル を，
＝ − ,すなわち， ＝ − として，
＝ ,
, ＝ ,0 0,−
( 0のとき）
( ＜0のとき）
により与えられる。もちろん，これは⑶の基底解 というだけでなく，主可能基底解である。なお，
⑸で定義した被約費用ベクトルと，ここで用いた 残差ベクトルとは互いに矛盾しない点に注意して いただきたい。また，基底 は 個の観察の説 明変数であり，対応する係数ベクトル は，超平 面 , ＝ が 個の点 , ∈ を 通るように定められる。
以上のように，LP ⑷の任意の基底 が与え られれば，それに対して，双対可能基底解のうち の非基底変数の値が容易に定められる。基底変数 の値は⑹により定められる。なお，例えば，
＝ を βの最小２乗推定量（OLS 推 定量）として，最小２乗残差ε＝ − の絶 対値 εを昇順に並べ，先頭から 個に対応する 観察の番号の集合を初期基底 とすることができ る。
4. ブートストラップ法による分散の推定と 初期可能基底
LADE の漸近分布は撹乱項の密度関数に依存 するため，推定量の分散を推定することが容易で はない。そのため，ブートストラップ法を用いる ことが効果的である（Greene［5, 16.3.2］）。
線形回帰モデルのためのブートストラップ法とし ては，
1. 説明変数と被説明変数のペアについてのブ ートストラップ（paired bootstrap）
2. 残差についてのブートストラップ（boot- strap based on residuals）
3. external bootstrap ま た は wild boot- strap
の３つがよく知られている（Shao and Tu ［7, 7.2.2］）。どのブートストラップ法を採用すべき かは，モデルの想定に依存する。例えば，１は無 作為標本であることのみを前提とした場合，２は 回帰モデルの撹乱項が独立で同一の分布に従う場
合，３は回帰モデルの撹乱項の分布が同一でない 場合に採用される。
以下では，３つの方法それぞれに対応した LP を整理し，より効率的な算法を開発する可能性に ついて検討する。
1. Paired bootstrap
説明変数と被説明変数のペアについてのブート ストラップでは，文字通り， 個のペア , ,
, , , , か ら，繰 り 返 し を 許 し て，大 きさ の標本を抽出する。表１の第１列は， ＝ 5の場合の３通りのブートストラップ標本の例で ある。繰り返しを許した標本抽出なので，１回だ け抽出されている観察もあれば，複数回抽出され る観察，あるいは１回も抽出されない観察もある。
１つのブートストラップ標本中での抽出回数を数 え上げたのが，表１の第２列である。
以上のように，１つのブートストラップ標本中 での抽出回数 , , , を定義することで，ブ ートストラップ繰り返し計算の各回に解くべき LP は，これまでと同様に(4)としてあらわされ る。制約式の右辺ベクトル と有界変数の上 限2 がブートストラップ標本に依存して変化す るので，元の標本に基づいて得られた最適基底解 は，一般に主可能ではなく，また，双対可能でも ない。したがって，ブートストラップ繰り返し計 算の各回毎に，前節で述べた方法により初期可能 基底解を求める必要がある。
2. Bootstrap based on residualsと external bootstrap
残差についてのブートストラップでは，元の標 本に基づいて得られた残差から再抽出を行う。元 の標本に基づいて得られた推定量を ，残差 をε＝ − ＝1, , と す る。こ の 残 差 ε,ε, ,εから繰り返しを許して抽出された大き さ の標本を ε,ε, ,ε とすると，ブートス 表１ Paired bootstrapの例：標本の大きさが
５の場合
ブートストラップ標本 (w,w,w,w,w) (2,1,4,2,1） (2,2,0,1,0）
(1,4,3,5,3） (1,0,2,1,1）
(3,5,3,2,4） (0,1,2,1,1）
ト ラ ッ プ 標 本 は , ＋ε , ,
＋ε , , , ＋ε である。すなわち，ブ ートストラップ繰り返し計算の各回に解くべき LP は，これまでと同様に⑷としてあらわされる。
元の標本とブートストラップ標本とを比較すると，
説明変数行列 とウェイト は変化せず，目的関 数ベクトル のみが変化する。この 性 質 は，
external bootstrap にも共通である。制約式と有 界変数の上限は元の標本の場合と等しいので，元 の標本に基づいて得られた最適基底解は，一般に 双対可能ではないが，主可能である。したがって，
ブートストラップ繰り返し計算の各回毎に２段階 単体法を適用する必要はなく，元の標本に基づい て得られた最適基底解を初期基底解として，主単 体法に基づく算法を適用することができる。
他方，双対単体法に基づく算法では，元の標本 に基づく LP の最適基底解が双対可能でないから，
前節末に述べた OLS 残差を用いる方法で可能基 底を求める必要がある。あるいは，基底だけは元 の標本に基づいて得られた最適基底を用いて，非 基底変数の値を変更することで可能基底解を求め ても良い。
以上の考察から，paired bootstrapの各ブート ストラップ標本に対しては，元の標本に適用する のと同じ算法の適用が必要であることが確認され た。他方，bootstrap based on residuals（およ び external bootstrap）では，元の標本に基づい て得られた最適基底解をブートストラップ繰り返 し計算において活用できることが確認された。し たがって，次節の計算実験では bootstrap based on residualsのみを考察の対象とする。
5. 計 算 実 験
付録にある通り，計算実験に必要なコードを MATLAB の M-script/functionsと し て 作 成 し た。単体法の計算効率はピボット演算の効率に大 きく依存するが，この点についての工夫は施して いない（むしろ，プログラムの見通しを優先する ために，非効率な演算を採用しているとさえ言え る）。そのため，以下では計算時間（CPU 時間）
はあくまでも参考値として示すものであり，繰り
返し回数とピボット演算回数を基準に各算法の計 算効率を比較検討する。
以下のデータ発生過程について計算を行う。
＝ , , , ＝1, , ,
＝1, 0,1 ＝1, , , ＝2, , , η ,ε＝η 3 ＝1, , ,
＝ 24
−1∑ ＋ε ＝1, , . 撹乱項は標準化された自由度３の 分布に従い，
定数項以外の説明変数は区間（0,1）上の一様分 布 に 従 う。ま た，理 論 的 な 決 定 係 数 は 2÷3＝
0.6666…である。擬似乱数は MATLAB のサブ ル ー チ ン randを 用 い る。標 本 の 大 き さ は ＝ 100, 200, 500の３通り，定数項を含む説明変数 の個数は ＝5, 10, 20の３通りについて，それ ぞれを 50回ずつ繰り返して計算実験を行った。
結果（計算時間，繰り返し回数，ピボット演算回 数 の 50回 の 平 均）は 表 ２ の 通 り で あ る。算 法
「主」は，主単体法に基づく有界変数法であり，
初期可能基底解が得られていないため，２段階法 を用いている。算法「双対」は，双対単体法に基 づく有界変数法であり，初期可能基底解は最小２ 乗残差の絶対値の小さい観察を選んで作成する。
表２ 計算実験の結果：主単体法と双対単体法に 基づく算法の比較
算法 n k 計算時間
〔秒〕
繰り返し 回数
ピボット 演算回数 主 100 5 0.221 153.1 111.7 双対 100 5 0.043 18.0 18.0 主 200 5 0.476 294.0 208.4 双対 200 5 0.062 22.7 22.7 主 500 5 1.528 713.0 499.0 双対 500 5 0.163 42.9 42.9 主 100 10 0.301 188.4 148.8 双対 100 10 0.085 35.1 35.1 主 200 10 0.666 361.0 279.2 双対 200 10 0.140 47.5 47.5 主 500 10 2.203 873.1 663.8 双対 500 10 0.309 73.4 73.4 主 100 20 0.444 233.8 193.4 双対 100 20 0.168 59.7 59.7 主 200 20 0.991 454.1 370.4 双対 200 20 0.268 82.0 82.0 主 500 20 3.407 1091.2 875.8 双対 500 20 0.783 160.7 160.7
いずれの算法も，標本の大きさ ，説明変数の個 数 の増加に伴って演算回数が増加する。算法
「主」では，非基底変数を０から上限（または，
その逆）に変更するだけでピボット演算を伴わな い繰り返しが行われ得る。繰り返し回数（ピボッ ト演算回数を内数として含む）は，ピボット演算 回数の 1.2倍から 1.4倍程度である。
標本の大きさ と説明変数の個数 のいずれの 組み合わせについても，算法「双対」の方が算法
「主」よりも効率が良い（ピボット演算回数が少 ない）。両者の効率が最も近い場合（ ＝100最小 値， ＝20最大値），「主」によるピボット演算回 数は，「双対」の約４倍である。最も差の大きい 場合（ ＝500最大値， ＝5最小値）では，この 比率は 18倍近い。
次に，ブートストラップのための算法について の比較検討を行う。計算実験の結果は表３の通り である。比較対象の３つの算法のうち，算法「双 対」は上で検討したものと同じである。算法「主
‑1」は，主単体法に基づく算法であるが，２段階 法ではなく，元の標本から得られた最適基底解を 初期可能基底解として用いる点が，算法「主」と の相違点である。算法「双対‑1」では，算法「主
‑1」と同じ基底を初期値として用いるが，双対可 能とするために算法「双対」と同様に非基底変数 の値を定める。
算法「主‑1」と「双対‑1」とを比較すると，両 者の効率が最も近い場合（ ＝100最小値， ＝ 20最大値），「主‑1」によるピボット演算回数は，
「双対‑1」の約 1.2倍であり，最も差の大きい場 合（ ＝500最大値， ＝5最小値）では，約 1.7 倍である。上で述べた「主」と「双対」との比較 と同様に，標本の大きさが大きいほど，説明変数 の個数が小さいほど，「主‑1」と「双対‑1」の効 率は乖離する。しかし，効率の相違はたかだか 1.7倍程度であり，「主」と「双対」の相違が最 大で 18倍であったことと比べると，驚くほど近 いと言える。参考値として載せた計算時間につい ては，説明変数の個数 の増加に伴って効率が逆 転するほどである。このことから，「主」が「双 対」と比べてかなり非効率であったのは，２段階 法により初期可能基底解を求める過程（フェーズ １）の効率の低いことが原因であったと思われる。
すなわち，元の標本についての最適基底解（元の 標本について LADE を求める際には利用できな いが，ブートストラップ繰り返しの各回には利用 できる情報）を活用して初期基底を与えることで，
主単体法に基づく有界変数法の計算効率が劇的に 改善されることが確認された。
次に，算法「双対‑1」と「双対」とを比較する。
両者の相違は初期（双対）可能基底解の相違のみ である。全体的に「双対」の方が効率が良く，効 率の比は２倍から８倍程度である。このことから，
期待値に相当する元の標本から得られた LADE よりも，各ブートストラップ標本から得られる OLSE の方が，計算効率向上に資する初期値で あることが分かる。
表３ 計算実験の結果：ブートストラップのための 算法の比較
算法 n k 計算時間
〔秒〕
繰り返し 回数
ピボット 演算回数 主‑1 100 5 0.125 74.5 72.8 双対‑1 100 5 0.116 48.3 48.3 双対 100 5 0.044 18.5 18.5 主‑1 200 5 0.270 145.1 141.8 双対‑1 200 5 0.259 93.5 93.5 双対 200 5 0.072 26.2 26.2 主‑1 500 5 0.869 347.4 339.4 双対‑1 500 5 0.837 216.9 216.9 双対 500 5 0.149 38.9 38.9 主‑1 100 10 0.169 94.1 93.9 双対‑1 100 10 0.177 70.7 70.7 双対 100 10 0.085 34.2 34.2 主‑1 200 10 0.390 184.3 183.9 双対‑1 200 10 0.411 137.4 137.4 双対 200 10 0.146 48.8 48.8 主‑1 500 10 1.284 442.6 441.3 双対‑1 500 10 1.376 324.4 324.4 双対 500 10 0.319 75.9 75.9 主‑1 100 20 0.244 116.1 116.1 双対‑1 100 20 0.276 98.6 98.6 双対 100 20 0.164 58.6 58.6 主‑1 200 20 0.588 241.5 241.5 双対‑1 200 20 0.663 202.0 202.0 双対 200 20 0.324 98.9 98.9 主‑1 500 20 2.095 588.7 588.6 双対‑1 500 20 2.354 482.3 482.3 双対 500 20 0.793 163.1 163.1
6. ま と め
計算実験により，最小絶対偏差推定量（LADE) を求めるための算法としては， 主単体法より も双対単体法の方が効率的であること， 主単 体法の非効率さは初期可能解の与え方（または，
求め方）に大きく依存すること， ブートスト ラップ繰り返しでは，元の標本に基づく最適基底 解を初期値として利用して，２段階単体法のフェ ーズ１を回避することで計算効率を大きく改善で きること，および 最小２乗残差に基づく双対 可能基底解は初期値として振る舞いの良いこと，
が分かった。この結果は，撹乱項が自由度３の 分布に従うものとして行った計算実験から得られ たものである。より裾の厚い分布，歪みのある分 布など，最小２乗推定量と LADE との乖離が大 きいと予想されるときには，本稿とは異なる振る 舞いを示すかもしれない。このような場合にも振 る舞いの良い初期値の与え方も含めて，LADE についてのブートストラップ法のための算法につ いて包括的な検討を行うことは，今後に残された 課題である。
［謝 辞］
本稿について，西郷浩氏（早稲田大学）から貴重なコ メントをいただいた。それらは改訂に際して大変有益で あった。記して感謝いたします。もちろん，本稿にあり 得べき誤りの責任は，筆者のみが負うものである。なお 本研究の遂行に際し，科学研究費補助金（奨励研究(A) 12730020）の助成を受けた。
［参 考 文 献］
［1］ 今野浩（1987），『線形計画法』日科技連。
［2］ Amemiya, T.(1985),Advanced Econometrics.
Cambridge, MA:Harvard University Press.
［3］ Armstrong, R. D., and M. T. Kung(1982), “A Dual Algorithm to Solve Linear Least Absolute Value Problems,” Journal of Operational Research Society33, 931‑936.
［4］ Charnes, A., W. W. Cooper, and R. Ferguson
（1955), “Optimal Estimation of Executive Com- pensation by Linear Programming,”Mathemati- cal Science1(2), 138‑151.
［5］ Greene, W. H.（2003),Econometric Analysis,
5th ed. Upper Saddle River, NJ:Prentice-Hall.
［6］ Maros, I.（2003), “A Generalized Dual Phasee-2 Simplex Algorithm,”European Journal of Opera- tional Research149(1), 1‑16.
［7］ Shao, J., and D. Tu（1995),The Jackknife and Bootstrap. Springer-Verlag, New York.
付録 A 計算実験に使用した M ATLAB プ ログラム
以下，計算実験に使用した MATLAB プログ ラム（印刷の都合上，行の折り返しなどについて 若干の編集を施したもの）を示す。
A.1 計算実験全体を行うための M-script
A.2 主単体法に基づく有界変数法のための M- function
A.3 双対単体法に基づく有界変数法のための M-function
