どのようにして分位点回帰の推定は線形計画法でなしうるのか

食費を節約したいし痩せたいけどついついコンビニでドリアとか買っちゃう社会の闇

Introduction

分位点回帰モデルに関する解説を読むと、「推定のところでは線形計画法を使う」と書いてある（Koenkerの有名な解説書等）。
分位点回帰の推定において、何が最小化されているかというのはそれなりにとっつきやすいんだが、
それがなぜ線形計画法で解けるのかということは何となくわかった気になってそのままにしていた。

なぜか急に気になって電車の中でスマホで調べたら、
以下のStackExchangeでの質問への回答がなかなかにわかりやすかった*1ので、自分用にメモしておく。

stats.stackexchange.com

線形制約問題自体については、最近たまたま見つけた以下の岩波講座「最適化法」の第2章（著者：今野浩＝ヒラノ教授）がこれまで読んだ中で一番過不足なく説明されていると感じた。さすが第一人者。

岩波講座応用数学〈15〉〔方法7〕最適化法／〔方法8〕離散最適化法とアルゴリズム

作者: 藤田宏,今野浩,田邉國士,茨木俊秀
出版社/メーカー: 岩波書店
発売日: 1998/08/03
メディア: 単行本
この商品を含むブログを見る

前提

線形計画法の標準形

線形計画問題の標準形は、以下のようにあらわされる。

$min \ \ c^T z , Az = b$

ここで、 $z ,\ c , \ b$ はいずれもk × 1、 $A$ はn × kの行列であるとする。
また、 $z >=0$ となっている。

解き方（単体法、内点法など）はともかく、この形で示すことができれば、
線形計画法として解にたどりつく（無限解や解が存在しないという判定も含め）ことができる。

分位点回帰における「最小にすべきもの」

一般に、回帰モデルとは、個人iに関する従属変数 $y_i$ を独立変数と回帰係数の線形結合と誤差の和として表すものである。
$y_i = \beta^T x_i + \epsilon_i$

さてここで、教科書的な説明を振り返っておくと、通常の線形回帰モデル（OLS）においては、
「誤差の二乗和を最小化する」という条件を満たすように回帰係数が求められる。すなわち、
$\beta \in R^k , ( s.t. \ \ min \ \Sigma_i^N \epsilon_i^2 )$

となる。一方、分位点回帰モデル（以下QRMと略記）においては、「誤差の重み付き和を最小化する」という条件を満たすように回帰係数が求められる。すなわち、

$\beta \in R^k , ( s.t. \ \ min \ \Sigma_i^N \rho_i |\epsilon_i |)$

ここで、weightである $\rho$ は以下のように定義される
$\rho_i = \tau \ ( \rm{if} \ \epsilon_i >=0 ) , 1 - \tau \ ( \rm{if} \ \epsilon_i < 0 )$

ここで $\tau$ は、独立変数が与える影響を推定したい分位点であり、 $0 < \tau < 1$ である。
よりわかりやすく日本語で説明すると、各ケース $i$ について「観測値 - 予測値」の符号によって重みを変えているということである。
たとえば20%点に対してのQRの場合、観測値が予測値を超える場合の重みは0.2、そうでないときは0.8となる。

何がわからなかったのか

これまで何がわからなかったのか、すなわち「分位点回帰は線形計画法で解ける」と言われてもどこが引っ掛かっていたのかというと、

そもそも回帰係数 $\beta$ が決まらないと $\epsilon_i$ の符号は決まらないのに、「 $\epsilon_i$ を最小化するように係数を推定する」というのは困難なのでは？（相互依存性の問題1）
$\epsilon_i$ の符号が決まらないとウェイト $\rho_i$ も決められないが、 $\beta$ が未決定だと前述のように $\epsilon_i$ もわからないのでウェイトがどうなるのかは事前にわからないのでは？（相互依存の問題2）

ということで、つまるところ
回帰係数 $\beta$ の推定値 ⇔ $\epsilon_i$ の値 ⇔ $\epsilon_i$ の符号　⇔　ウェイト $\rho_i$ の値
という相互依存関係をどのように扱うのかということがわからなかったということである。

相互依存性への解決法としての同時的決定

今回勉強してみての結論からいうと、上記の問題は
「 $\epsilon_i$ 、 $\beta$ 、 $\rho_i$ を同時決定されるパラメータとして解く」
によって解消される。

線形制約側の話

まず線形制約側について考える。回帰モデル
$Y_i = \beta^T X_i + \epsilon_i$
をさきほどみた線形制約 $Az=b, z >=0$ の形に変えることを考える。
とりあえず縦にn個つなげて、
$Y = \beta^T X + \epsilon$
とする、左辺右辺ともに要素n個のベクトルとなる。

さてここで行列Aにあたるのは、観測された独立変数 $X$ となる。
「変数」という言葉で惑わされそうになるが、観測値は固定なので、スカラの制約行列となるのである。
そして線形制約式の右辺 $b$ にあたるのは、従属変数列 $Y$ である。

したがって、線形制約式における変数 $z$ にあたるのは、残った二つ。
係数 $\beta$ と、誤差 $\epsilon$ である。

しかし、これらは定義上負の値をとりうる（これを線形計画法の用語で「自由変数」という）ので、
以下のように二つの非負変数ベクトルの差として表現しなおすこととする。
$\epsilon = u - v , \ (u >=0 ,\ v >= 0 )$
$\beta = \beta' - \beta'' ( \beta' >=0 , \ \beta'' >=0 )$

すると、
$Y = \beta^T X + \epsilon$
$\iff \ Y = ( \beta' - \beta'' )^T X + u- v$
であり、ここで次元nの単位行列 $I$ を導入すると
$\iff \ Y = ( \beta' - \beta'' )^T X + (u- v)^T I$
となるから結局、
$A = \begin{pmatrix} X & -X & I & -I \end{pmatrix}$
$z = \begin{pmatrix} \beta' & \beta'' & u & v \end{pmatrix}$
とおくと、
$Az = Y$ という綺麗な線形制約の形になっているのである。

最小化対象式側の話

残りは $\textrm{min} \ c^Tz$ のほうだが、
変数zは先ほど定義したので、c（スカラのベクトル）を適切に定義すればいい
$c = \begin{pmatrix} \textbf{0} & \textbf{0} & \tau 1_N & (1- \tau) 1_N \end{pmatrix}$
ただしここで $1_N$ はN×1ですべての要素が1の行列である。

cの定義から明らかであるが、最小化対象（＝cとzの内積）には、回帰係数が出てこないようになっている。
実測値＞予測値となるときは $\epsilon_i = u_i$ 、実測値＜予測値となるときは $\epsilon_i = - v_i$ となる*2ので、
$c^T z = \begin{pmatrix} \tau 1_N & (1- \tau) 1_N \end{pmatrix}^{T}$ $\begin{pmatrix} u & v \end{pmatrix}$
は最小すべき重み付き絶対値和 $\Sigma_i^N \rho_i |\epsilon_i |$ になっている。

あとは、これを普通の線形計画問題として解けばQRの推定値が得られるのである。
めでたしめでたし。

最後に

なんか「一見相互依存に見えていて（デッドロックを起こしてしまい）解けない問題を同時決定問題として解く」というのは、
統計学や工学など理系の問題だけでなく、社会科学や、政策決定などのより一般的な問題解決においても適宜必要な知恵ではないかと思った。
囚人のジレンマがジレンマなのは同時決定できる超越的な主体が存在しないからだよね、って話である。

線形計画もわかったような気になってまた忘れるを繰り返している気がするが
時間があったら前述の岩波講座の忘れちゃいけない部分だけまとめときたい。

Enjoy!

*1:なのに質問者は「長くてわからん」って言っててああ無情。

*2:厳密にはこれは別途証明する必要があるが省略する