科学計算・機械学習ソフトウェアのドメイン駆動設計

目次:

はじめに

　この記事は、ドメイン駆動設計に親しみのあるソフトウェアエンジニアで、科学計算・機械学習分野のソフトウェア設計に興味のある人に向けて書いています。機械学習、統計解析やそのためのデータ処理の経験がないと、これらの分野特有のメンタルモデルを認識するのにちょっとしたハードルがあると思ったので、そうした認識のギャップを埋めるのに役立てばと思い記事を書くことにしました。特に NumPy や PyTorch, Polars といったライブラリに触れてプロダクトコードの設計にどう組み込もうかと考えている方への一助となればと思います。

ドメイン駆動設計とは

　ドメイン駆動設計にある程度馴染みのある方に向けて書いているので、DDD自体の説明は簡単な (筆者の拙い理解による) ものに留めます。質の高い説明はEric Evansの書籍などをご参照頂ければと思います。非常にざっくり言えば、DDDとは 「ソフトウェアの振る舞いを、技術的詳細に依存しないビジネスロジックとみなしてドメインオブジェクトとその操作で表現する」 ことであり、 「エンジニアがビジネスメンバーと一体になってそうしたドメインモデルを深化させていくことをエンジニアリング活動の中心に据える」 ことであると筆者は理解しています。プロダクト企業の開発組織で働いたことのある筆者の感覚としても、DDDを実践していくことでプロダクトが成長していく方向性とコードベースが成長していく方向性を揃えていくことが出来て、開発生産性を高く保ちビジネスの持続的成長にコミットできるという強い実感があります。

コレクションを操作するときのメンタルモデル

　多くのソフトウェアエンジニアが目にするであろう、DDDで書かれたKotlinのコード例を書いてみます:

data class Product(val name: String, val price: Double) {
    fun applyDiscount(rate: Double): Product {
        // 価格に割引率を適用した新しいProductを返す
        return this.copy(price = this.price * (1 - rate))
    }
}

class Products(private val items: List<Product>) {
    fun applyDiscountToAll(rate: Double): Products {
        // コレクション内の各Productに割引を適用
        val discountedItems = items.map { it.applyDiscount(rate) }
        return Products(discountedItems)
    }
}

簡略化のために name, price や rate などのドメインをプリミティブ型とした点はご容赦ください。 Product というドメインがあり、その first class collection としての Products というドメインを記述しています。

　ここで注目してほしいのは applyDiscountToAll メソッドです。 map メソッドを使って、「ひとつひとつのプロダクト全てに割引を適用する」という表現になっています。これが「全てのプロダクトに割引を適用する」というビジネスロジックに対して普通の人が抱くメンタルモデルではないでしょうか。テスト容易性が高く、単一責任の原則が守られ、自己説明的なコードになっていると言えるでしょう。

　他方、特に大量データ処理を要する科学計算や機械学習分野においては、ちょっと違った考え方をする必要があります。その考え方の理解に必要な背景知識として、次セクションではSIMDと呼ばれるCPUの機能について紹介します。

SIMD (Single Instruction, Multipe Data)

SIMD (ベクトル化演算) とは

　SIMDとは、読んで字の如く、一度の命令で複数のデータを一度に処理出来るというCPUの機能のことです。配列(ベクトル)データに対する同一の演算を一度に処理出来ることから ベクトル化演算 とも呼ばれ、機械学習分野ではこちらの呼び名の方が馴染みがあると思います。昨今のLLMやAIの文脈で出てくる「ベクトル」とは全く別の概念なので注意してください。

　CPUにはベクトルレジスタと呼ばれる小さなメモリ領域がついています。例えば256bitのベクトルレジスタ幅であれば、64bit整数を256 ÷ 64 = 4つ格納することが出来ます。同じ演算を4つのデータに同時に適用出来るので、SIMDを活用出来ればこの場合は理論上4倍計算が早くなります。SIMDは厳密には並列処理の一種ですが、マルチコアで並列処理をしているのではなく、あくまでシングルコアの動作として複数データを同時に処理出来る、というわけです。みなさんのデスクトップPCや各社クラウドサービスが提供するサーバーには、大抵の場合256bit以上のベクトルレジスタ幅を持つCPUが搭載されています。

ベクトル化演算を使った実装

　試しに、上記の割引を適用する実装を Python の NumPy パッケージを使って実装してみます:

import numpy as np

prices = np.array([100, 90, 120, 100])
rate = 0.2
discounted_prices = prices * (1 - rate)

簡略化のためにプロダクトではなく価格のみのコードにしました。一見すると、Kotlinのコードと全く同じ map 処理をしているように思えます。しかし、NumPy を使った実装により、シングルコアで実行されているにもかかわらず割引の計算が配列の各要素に対して同時に実行されている、という点がKotlinの実装とは異なっています。

　なぜ上記のNumPyのコード例ではベクトル化演算を実行出来ているのでしょうか。prices オブジェクトは、numpy.ndarray というNumPyが提供する配列の型を持っています。この ndarray の要素は Python のプリミティブ型としての整数ではなく、これまたNumPyが提供する整数の型を持っています:

>>> type(100)
<class 'int'>

>>> type(prices.dtype)
<class 'numpy.dtypes.Int64DType'>

さらに、* や - などの演算子は ndarray クラスに実装されたメソッドのシンタックスシュガーです。つまり配列に対する四則演算の処理はNumPyに実装されているのです。NumPyそのものはC言語で書かれていて、その中に配列に対するベクトル化演算が実装されているために、先程のコードでベクトル化演算が可能だったというわけです。四則演算に限らず、配列の変換処理やより複雑な計算もNumPyがAPIとして提供しているものであればベクトル化されています。

科学計算のドメインはフレームワークに依存する

ベクトル化演算のメンタルモデル

　機械学習やそのためのデータ処理ではパフォーマンスが求められるため、プログラミング言語が提供するプリミティブなループ処理でデータを操作することは無く、上述のフレームワークがサポートするベクトル化演算をフル活用して実装されます。データサイエンティストや機械学習エンジニアは、ベクトル化されていない科学計算の実装に対して強い "code smell" を感じるように訓練されています。上記のKotlinコードの例のように「ひとつひとつのデータを順々に処理する」という感覚は無く、「データのカタマリを同時に処理する」 というハードウェアの動作が透けて見えるような目線を持ってデータ処理を捉えているのです。

　NumPyで実装されていたベクトル化演算の状況を改めて整理してみましょう:

• データの入れ物 (配列) やその要素はNumPyの型
• データを操作するロジック (API) はNumPyの実装

データの入れ物もロジックもフレームワーク依存なのです。そのためドメイン設計をするにあたっても、ドメインオブジェクトやその操作がフレームワークに依存することになります。

フレームワークの設計思想

　データサイエンティスト・機械学習エンジニアは科学計算 (機械学習) のライブラリやデータ処理用のフレームワークを頻繁に用います:

• 科学計算: NumPy, PyTorch, LightGBM
• データ処理: Pandas, Polars

実はこれらにはあるひとつの共通した設計思想があります。いずれも「フレームワークが提供するデータ型とAPIの世界の中で、ベクトル化演算を実現する」という作りになっているのです。特にPandasやPolarsなどは行ごとのデータ処理を行わずに、列を単位としたベクトル化演算がサポートされているため、 列志向 (columnar, column oriented) データ処理 のフレームワークと呼ばれています。

フレームワークにまつわる誤解

　ソフトウェアエンジニアがフレームワークと聞くと、webフレームワークのように「クリーンアーキテクチャの同心円の最も外側に位置すべき技術的詳細」というイメージを抱くと思います。ですが科学計算 (や機械学習) 分野ではフレームワークが担う役割はwebのそれと異なります。データの入れ物やロジックの実体を担っているため、フレームワークを操作することがビジネスロジックの中心になるのです。フレームワークはインフラ層にあるべき技術的詳細ではなく、プリミティブな型とAPIを提供するプログラミング言語そのものに近い性格を持っているとも言えるかもしれません。

　この「ドメインやドメインを操作するロジックがフレームワークに依存している」というのは私独自の考えであり、多くの機械学習エンジニアが同意する標準的な意見ではないと思います。標準的な考え方としてはどちらかといえば、データ処理バッチや推論APIなどの計算負荷の高いコンポーネントを個別のマイクロサービスに切り出して、その中ではドメイン設計を明示的に考えず手続き的な実装とする、というのが主流な気がします。ですがフレームワークの位置づけがwebとは異なるという点は間違いないので、知っておいて損はないと思います。

ドメインオブジェクトの実装例

　ごく簡単にですが、このような考え方に則ったドメインオブジェクトの実装例を載せてみます。ある程度計算負荷のある集計処理を要するデータサイズの大きいドメインを扱う例として、金融のドメインにおける約定 (execution) のデータを扱う例を考えてみました:

import polars as pl
from pydantic import BaseModel, ConfigDict, field_validator

EXECUTION_SCHEMA = {
    "side": pl.Enum(["BUY", "SELL"]),
    "amount": pl.Float64(),
    "symbol": pl.Enum(["USD/JPY", "EUR/JPY"]),
    "book_id": pl.Enum(["1", "2", "3"]),
}


class Execution(BaseModel):
    data: pl.DataFrame

    model_config = ConfigDict(arbitrary_types_allowed=True)

    @field_validator("data", mode="after")
    @classmethod
    def validate_schema(cls, v: pl.DataFrame) -> pl.DataFrame:
        expected_schema = pl.Schema(EXECUTION_SCHEMA)
        if sorted(expected_schema.items()) != sorted(v.schema.items()):
            raise ValueError(
                f"Schema mismatch. Expected: {expected_schema}, Got: {v.schema}"
            )
        return v

    def is_empty(self) -> bool:
        return self.data.height == 0

　PythonでDDDを実践する際は、Pydanticで型ヒントによるバリデーションを導入するのが現代的で標準的な実装です。約定データの実体は、データ処理ライブラリのデファクトである Polars のデータフレームオブジェクトとなっています。データフレームを Execution というドメインオブジェクトでラップすることで、例えばデータの実体をPandasに移行したとしてもユースケース(アプリケーション)レイヤーから見るとソフトウェアの挙動は不変にすることが出来ます。上記では is_empty しかビジネスロジックがありませんが、 Execution に何かロジックを持たせたいときは、別のドメインオブジェクトを返すメソッドを実装してあげれば良いでしょう。

　データフレームのスキーマを表すオブジェクトとして EXECUTION_SCHEMA を導入しています。Execution の初期化時に validate_schema メソッドが実行され、 EXECUTION_SCHEMA に基づいたデータフレームの型チェックがランタイムで行われる仕組みです。バリデーションによるオーバーヘッドが気にならなければ、データフレームのバリデーションツールである Pandera を導入するとよりシンプルで可読性の高い実装が出来ると思います。

チームとしてドメインを育てる

　このようにフレームワークに依存したドメイン設計は、DDDの基本的な考え方からは逸脱しています。アジャイル開発における守破離でいう「破」に該当する考え方だと思います。機械学習などの科学計算におけるドメイン設計には、「こう考えるのが標準的だ」という守にあたるような正解はありません。そのため、チームの中で自分たちなりの正解を探す試行錯誤を続ける気概を持ち、「私たちはこのドメインをこう表現しよう」という議論を重ねて合意形成していかないと、機械学習システムを通じてビジネスの要求に応え続けることは出来ないと思います。私が示した考え方もあくまで選択肢のひとつであり、もっと良い正解を探していきたいという前提で書いています。本記事を読んだソフトウェアエンジニアにデータサイエンティストや機械学習エンジニアの視点を知ってもらえて、機械学習システムの改善点を考えるきっかけになってもらえたら幸いです。

ラベル割り当ての層化傾向スコアを用いた共変量シフト下の selective prediction

目次:

はじめに

機械学習システムの継続的デリバリ

　機械学習システムの開発において、継続的デリバリの実現は重要な課題の一つです。推論結果の継続的な品質向上という観点では、推論対象全てのインスタンスについて推論結果をリリースし、継続的デリバリにより逐次的に推論対象全体の精度を改善していく、というリリース戦略が最も典型的でしょう。すなわち、推論結果のカバレッジが100%であることを保証しながら精度を逐次的に上げていくというアプローチです。

　他方、「カバレッジは低くても良いから、一定の精度以上の推論対象のみリリースしたい」というケースもありえます。医療ドメインや自動運転などミッションクリティカルな応用が特にこれに当たります。このような、予測モデルが予測結果に自信を持てないときに予測結果を返さない "reject option" がありえる予測のことを、 "selective prediction" と呼びます*1。この場合、推論対象のサブセットの精度が一定以上であることを保証しながら、そのような推論対象のカバレッジを広げていくアプローチになります。

図 (1): 推論対象のカバレッジと精度の逐次的な改善。筆者作成。

　図 (1) は上記のようなリリース戦略の違いを表しています。あるモデルが実現可能な推論対象のカバレッジと精度の組み合わせを等高線として描画し、この等高線を右上方向へシフトさせていくことがモデルの改善であると考えます。通常は左図のように推論対象全体について精度を上げていくリリース戦略が採られますが、一定以上の精度となる推論対象のサブセットをリリースし、そのカバレッジを継続的に向上させるリリース戦略を示しているのが右図です。本記事では後者のようなカバレッジを上げていくスタイルの継続的デリバリに役立つ selective prediction の手法について考えます。

　ここでは、Selective prediction の問題設定を簡単に紹介しておきます (Geifman & El-Yaniv, 2019)。 を特徴空間、 をラベル空間とし、 を における分布とします。ラベルを予測する "prediction funciton" を とすると、 の に関する真のリスクは、 ] と表すことができます。また、 を reject option を使うかどうかの "selection function" とします。この２つを組み合わせた "selective prediction function" は、以下のように表せます:

\begin{align} (f,g)(x) := \begin{cases} f(x) & \text{if } g(x) = 1\\ \text{null} & \text{if } g(x) = 0. \end{cases} \end{align}

このような selective risk は以下のように定義できます:

\begin{align} R(f,g) := \dfrac{ E_{P(X,Y)} [\ell (f(x),y)g(x)] }{E_{P(X,Y)}[g(x)] } \end{align}

Selective risk が一定の閾値 以下という制約の下での最適な selective prediction function は、

\begin{align} \max_{\theta \in \Theta_h, h \in \mathbb{H}} & E_{P(X,Y)}[g_{\theta, h}(x)] \\ \text{subject to } & R(f_{\theta, h},g_{\theta, h}) \leq r. \end{align}

ここで と はそれぞれ structural risk minimization における仮説集合です*2。図 (1) の右図のプロセスは、学習やモデル選択により得られた をモデル改善を通じて最適な に近づけていくことに相当します。

ラベル割り当てと共変量シフト

　上記のように推論対象から一部のサブセットを選択することを考える上で重要なのが、共変量シフトの観点です。正解ラベルの割り当てメカニズムは通常は完全にランダム化されていることなく、ラベル付きデータと推論対象とで特徴量の分布が異なりうるため、ラベル付きデータで得られた精度が推論対象では成り立たないという問題が様々な機械学習ドメインの実務で頻繁に起こります。そのためラベル付きデータと特徴量の分布が類似する推論対象を選択する out-of-distribution (OOD) detection と呼ばれる技術が重要になりますが、特にテーブルデータの機械学習ドメインで、質的変数や欠損値を含むような高次元の特徴量の場合、分布の類似度を測るのが難しくなります。

ラベル割り当ての層化傾向スコアを用いた selective prediction

　共変量シフトが発生している (テーブルデータの) 機械学習タスクにおいて、特定の閾値以上の精度を満たす selective prediction モデルはどのように学習出来るでしょうか？本記事では、筆者が実務において考案した、ラベル割り当ての層化傾向スコアを用いる手法を紹介します。下記セクションよりこの手法を詳述しますが、コード例はこちらのノートブックを参照してください。

Adversarial validation

Adversarial validation とは

　本記事の手法は、adversarial validation と呼ばれる、共変量シフト下でラベル付きデータへの過学習を防ぐためのモデル選択の手法を参考にした点が多くあるため、はじめに adversarial validation について簡単に説明しておきます。

　各インスタンスに正解ラベルが付与されていれば 1、正解ラベルが付与されていない推論対象であれば 0 となる二値の変数を付与し、元の機械学習タスクで用いている特徴量でこの二値変数を予測する分類器 (adversarial classifier) を学習します。Adversarial classifier の精度が ROC AUC = 0.5 など random guess に近い状況であれば、ラベル付きデータとラベル無しデータの特徴量の分布はほぼ同一であると言えるでしょう。逆に adversarial classifier の精度が高い状況であれば、ラベル付きデータとラベル無しデータの特徴量の分布に差が生じてしまっていると言えます。このような adversarial classifier をモデル選択に利用する手法を adversarial validation と呼びます。具体的な利用の仕方は主に2通りあります。

　第一には、ラベル無しデータである推論対象に類似したラベル付きデータを選択し、これを精度評価に用いる手法です。シンプルにラベル付きデータである確率が特定の閾値以下のデータを精度評価に用いたり、傾向スコアマッチングや Inverse probability weighting を用いて推論対象に類似したラベル付きデータを評価データに選択するやり方などがあります (Pan et al., 2020)。

　第二には、特徴量選択に用いる手法です。adversarial classifier の特徴量重要度の高い特徴量を、元の機械学習タスクの特徴量から除くような特徴量選択をすることで、元の機械学習タスクにおいてラベル付きデータへの過学習を防ぐことが出来ます。

　Adversarial classifier は、質的変数や欠損値を含む特徴量であっても LightGBM や XGBoost などの GBM の手法を用いて学習出来るため、特にテーブルデータの機械学習において優れた手法と言えます。そのためKaggle などの機械学習コンペで古くから用いられてきました。

Adversarial validation の統計的因果推論的解釈

　あるインスタンスにラベルが割り当てられる確率を予測する adversarial classifier は、潜在アウトカムアプローチの統計的因果推論における傾向スコアの予測モデルであると解釈することが出来ます*3。

　インスタンス における ラベルの割り当てを 、特徴量を 、 ラベル割り当ての (真の) 傾向スコアを とします。Imbens & Rubin (2015) の 12.2.1 の "regular assignment mechanism" の仮定を置くと*4、Lemma 12.1 より、傾向スコアを条件づけるとラベル割り当ては特徴量と独立となります:

\begin{align} \newcommand{\ind}{\perp\!\!\!\!\perp} W_i \ind X_i | e(X_i). \end{align}

したがって、傾向スコアを条件付ければ、ラベルありデータとラベル無しデータで特徴量の分布は等しくなり、故にラベルありデータとラベル無しデータとで予測誤差の分布も等しくなります。直感的には、似たような傾向スコアのインスタンス群であれば、ラベルありデータで評価した元の機械学習タスクの精度がラベル無しのデータでも成り立っている、ということです。

　次に、adversarial classifier によりこのような傾向スコアを予測することを考えます。予測された傾向スコア は、式 (5) のようにラベル付きデータとラベル無しデータとで特徴量の分布がバランスするように予測されているはずです:

\begin{align} \newcommand{\ind}{\perp\!\!\!\!\perp} W_i \ind X_i | \hat{e}(X_i). \end{align}

　式 (5) の条件付き独立性が十分性を以て満たされているかを直接的に検定 (test) することは出来ませんが、必要性を検証 (assess) する実務的手段は確立されています。その代表的な手法のひとつが Imbens & Rubin (2015) や Lee (2011) でも紹介されているような傾向スコアの予測値の stratification (層化) と balancing test です。傾向スコアの予測値を 0〜0.1, 0.1〜0.2, ...といった "stratum" に分割し、stratum 内で特徴量がバランスしているかを記述統計量により検証するのです。「傾向スコアの予測値が balancing score として機能しているならば、少なくとも特徴量の記述統計量はラベル割り当てに依存していないはずだ」といったイメージです。

本手法のアルゴリズム

概要

図 (2): Selective prediction のアルゴリズム概要。筆者作成。

　図 (2) を用いてアルゴリズムの概要を説明します。目的変数を予測する学習モデルだけでなく、ラベルの割り当て有無を予測するモデル (adversarial classifier) を学習します。Adversarial classifier の out-of-fold prediction により得られたラベル割り当ての傾向スコアの予測値について、層ごとに balancing test を行います。さらに元の機械学習タスクについて目的変数を学習したモデルの精度評価を、傾向スコア予測値の層ごとに行います。Balancing test を通過した層と、精度の値が閾値以上となる層との共通集合が selective prediction の対象 (リリースの対象) となる、というアルゴリズムです。層化傾向スコアの balancing test が OOD detector として働き、層ごとの精度評価が selective predictor として働くようなイメージです。

層化傾向スコアの balancing test

アルゴリズム概要

　本セクションでは、前セクションで触れた balancing test として典型的に用いられる手法を紹介します。アルゴリズムはとてもシンプルで、各stratumで、各特徴量について、ラベル付きデータとラベル無しデータとでの特徴量単変量の平均値の差の検定を行います*5。Balancing test のアルゴリズムの疑似コードは以下のようなイメージです:

results = []
for stratum in strata:
    results_by_stratum = []
    for feature in features:
        result = mean_difference_test(labeled_data[feature], unlabeled_data[feature], significance=0.05)
        results_by_stratum.append(result)
    results.append(results_by_stratum)

実務上の考慮

　全ての特徴量について検定を行うと、検定結果に偶然の過誤が生じる多重検定の問題が発生するので、検定対象となる特徴量はなんらかのルールで絞ることが balancing test において一般的です。今回のユースケースでは、元の機械学習タスクと傾向スコアの双方に大きな影響を与えている特徴量がバランシングしていることが重要なので、双方の特徴量重要度の高い特徴量上位K個に対象を絞ると良いかと思います。

　因果推論の文脈ではほとんど言及がありませんが、テーブルデータの機械学習の場合は質的変数や欠損値を持つ特徴量も含む場合があります。欠損値を含む量的変数の検定には母比率の差の検定を用いて欠損の無い部分同士で平均値の差の検定を行ったり、質的変数の検定には Cochran's Q test (Campbell, 2021) を用いる方法があります。

層ごとの予測精度

　前セクションでは、層化傾向スコアを用いて各 stratum で特徴量の分布が類似しているかの balancing test を行いました。この stratum ごとにラベルありデータで予測精度を評価した場合、傾向スコアと予測精度はどのような関係になるでしょうか。傾向スコアが低いほど、ラベルありデータらしさが失われてラベル無しの推論対象のデータに近づいていくため、精度は下がっていくはずです。傾向スコアが高いほど精度は高くなりますが、傾向スコアが高くなりすぎるとラベルありデータの中でも外挿的なインスタンスが多くなってくるため逆に精度は下がっていくでしょう。したがって、精度の閾値を満たす strata は、傾向スコアの中央値 (およそ 0.5) よりもやや大きい方向にシフトして分布しているはずです。

Selective prediction

　上記の通り、傾向スコアの stratum には精度が高くなる領域が存在しているはずです。さらに、前セクションの balancing test が通っている stratum であれば、ラベルあり・無しデータで特徴量の分布が類似しているため、ラベルありデータで測った精度がラベル無しデータでもほぼ同じ値になっているはずです。すなわち、傾向スコアで層化した stratum ごとの精度のうちリリース基準を超えている stratum と、balancing test が通った stratum とがオーバーラップしている領域をリリース対象と出来るはずです。

図 (3): 層化傾向スコアの balancing test と評価指標の関係の概念図。筆者作成。

　図 (3) は、このような selective prediction における balancing test と精度の関係を表した図です。"balancing region" が balancing test を通過している stratum の範囲を表しており、"eval metric (above threshold)" がラベルありデータで評価したリリース基準を通過した stratum を表しています。両者がオーバーラップしている stratum に属する推論対象を selective prediction の対象とすることが出来ます。ラベル無しデータの中でもラベルありデータに類似したインスタンスで、かつ精度の閾値を満たすインスタンスのみが予測の対象となるようなイメージです。

特徴量選択によるカバレッジ最大化

　前セクションで、ラベル割り当ての層化傾向スコアを用いて selective prediction function を学習するアルゴリズムを紹介しました。本セクションでは、式 (3) に表される所与の selective risk の許容度の下でのカバレッジを最大化させるため、adversarial validation による特徴量選択を応用する方法を紹介します。

　カバレッジを向上させるひとつの方法は、ラベルつきデータとラベル無しデータの分布を近づけさせ、共変量シフトの影響を緩和させることです。Adversarial classifier の特徴量重要度などを用いて、傾向スコアの予測に寄与している特徴量を除くような特徴量選択を行います。これにより adversarial classifier の予測精度は random guess に近づいていき、balancing test は通りやすくなるはずです。

　一方、上記のように特徴量を除き続けていくと、元の機械学習タスクの精度が下がっていくことになります。すると精度 (またはselective risk) が閾値を満たす strata の数も少なくなっていき、逆にカバレッジは小さくなっていきます。

図 (4): Risk-coverage trade-off の概念図。筆者作成。

　最適なカバレッジは上記のトレードオフが釣り合うところに存在しています。 図 (4) は、このような risk-coverage trade-off と呼ばれる関係を示した図です。左図はラベルありデータでの精度は高いが balancing region の範囲が狭い状況を表しています。他方、右図は balancing region の範囲は広いが、ラベルありデータでの精度は低い状況を表しています。Adversarial validation による特徴量選択でこのトレードオフのバランスをとることができ、精度の閾値を満たしかつ balancing region となっているようなカバレッジが最も広くなるようにモデル選択をすることができます。

レファレンス

• Campbell, M. J. (2021). Statistics at square one. John Wiley & Sons.
• El-Yaniv, R. (2010). On the Foundations of Noise-free Selective Classification. Journal of Machine Learning Research, 11(5).
• Geifman, Y., & El-Yaniv, R. (2017). Selective classification for deep neural networks. Advances in neural information processing systems, 30.
• Geifman, Y., & El-Yaniv, R. (2019). SelectiveNet: A deep neural network with an integrated reject option. In International conference on machine learning (pp. 2151-2159). PMLR.
• Imbens, G. W., & Rubin, D. B. (2015). Causal inference in statistics, social, and biomedical sciences. Cambridge university press.
• Lee, W. S. (2011). Propensity score matching and variations on the balancing test. Empirical economics, 44, 47-80.
• Pan, J., Pham, V., Dorairaj, M., Chen, H., & Lee, J. Y. (2020). Adversarial validation approach to concept drift problem in user targeting automation systems at uber. arXiv preprint arXiv:2004.03045.
• Shalev-Shwartz, S., & Ben-David, S. (2014). Understanding machine learning: From theory to algorithms. Cambridge university press.

*1:Selective prediction は特に分類問題の文脈で研究が進んでいます (El-Yaniv, 2010; Geifman & El-Yaniv, 2017)。

*2: がハイパーパラメータなどのモデル選択における選択肢の集合であり、 が所与のモデルの下で学習可能なパラメータなどの集合を意味しています(Shalev-Shwartz & Ben-David, 2014)。すなわちここではモデル選択も含めた定式化をしています。

*3:この場合、ユーザーなどのアノテーターが自身の選好に基づいて最適化したい利益や効用が、統計的因果推論における潜在アウトカムに該当すると解釈することが出来ます。例えば映画の推薦タスクにおいてユーザーが自身の好みに基づいて映画のレーティングを行うとき、ユーザーによるレーティングが処置で、それによる心理的・経済的な便益がアウトカムであるとみなせるでしょう。潜在アウトカムをベースとした因果推論に興味のある読者は Imbens & Rubin (2015) などを参照してください。

*4:直感的には、「特徴量が同じ値ならばラベルの割り当てはランダム化されている」という仮定です。例えば、年齢という特徴量を使って身長を予測する機械学習タスクを考えます。このとき、10代以下の被験者にだけ正解ラベルをつけていて、推論対象が20代以上であるとします。すると傾向スコアは「10代以下ならば1、そうでなければ0」という deterministic な関数となってしまい、使い物にならなくなってしまいます。

*5:平均値の差の検定には、等分散性を仮定しない Welch の t 検定を用います

機械学習による区間予測入門 ⑥: 実装例

　本記事は機械学習による区間予測入門シリーズの6記事目です。他記事へのリンクはこちら: 

機械学習による区間予測入門 ①: はじめに

機械学習による区間予測入門 ②: 分位点回帰

機械学習による区間予測入門 ③: Conformal Prediction

機械学習による区間予測入門 ④: Conformalized Quantile Regression

機械学習による区間予測入門 ⑤: 予測区間の評価指標

目次:

はじめに

　区間予測の一連の記事を通じて、区間予測の代表的な手法や評価指標について紹介してきました。最後となる本記事では実際のデータを用いた実装例を記載します。まずはベンチマークとして LightGBM を用いた分位点回帰により精度を評価します。次に、分位点回帰で求めた区間を補正する手法として Conformalized Quantile Regression を実装します。CQR の実装では既存のフレームワークを用いず、筆者の設計の考え方も併せて紹介しました。

　まずは今回必要な依存関係をインポートしておきます:

from abc import ABC, abstractmethod
from typing import Any

import lightgbm as lgb
import numpy as np
import numpy.typing as npt
import plotly.express as px
import polars as pl
from pydantic import BaseModel
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import KFold
from tqdm.notebook import tqdm

データと評価指標

データ

　今回は scikit-learn を通じてインポート可能な California Housing dataset を用います。目的変数は地域ごとの平均住宅価格 (10万ドル単位) です。

housing = fetch_california_housing()
X, y = housing.data, housing.target

評価指標

　過去の予測区間の評価指標の記事に記載の通り、予測区間には、他の機械学習タスクと同様に正解データと予測値の双方を用いて算出する指標と、予測区間のみを用いて評価可能な指標の2種類があります。そのため2つのインターフェースを用意しました。

type NumericNDArray = npt.NDArray[np.number[Any]]


class PIMetricWithTarget(ABC):
    @abstractmethod
    def score(
        self,
        y_true: NumericNDArray,
        y_pred: NumericNDArray,
    ) -> NumericNDArray:
        raise NotImplementedError


class PIMetricWithoutTarget(ABC):
    @abstractmethod
    def score(
        self,
        y_pred: NumericNDArray,
    ) -> NumericNDArray:
        raise NotImplementedError

併せて、学習や予測、評価指標の算出時に必要なデータのバリデーターを用意しました。目的変数は1次元配列で、予測区間は上限と下限の2つから成る2次元配列となることをチェック出来るようにするためです。また、予測区間の下限が上限を超えないか否かもチェックします。

class BaseNDArrayValidator(ABC):
    @abstractmethod
    def __call__(self, a: NumericNDArray) -> None:
        raise NotImplementedError


class TargetVariableValidator(BaseNDArrayValidator):
    def __call__(self, a: NumericNDArray) -> None:
        if a.ndim != 1:
            raise ValueError("y_true must be 1-dimensional.")


class PredictionIntervalValidator(BaseNDArrayValidator):
    def __call__(self, a: NumericNDArray) -> None:
        if a.ndim != 2 or a.shape[1] != 2:
            raise ValueError(
                "y_pred[:,0] and y_pred[:,1] must be lower and upper bound of the prediction regions, respectively."
            )
        if not np.all(a[:, 0] <= a[:, 1]):
            raise ValueError("y_pred[:,0] must be less than or equal to y_pred[:,1].")

今回評価に用いる指標は、区間の精度を示すPICPと、区間の(正規化された)広さを示すNMPIWの2つとします。PICPには正解データが必要となりますが、NMPIWは予測区間のみから算出出来ます。指標の詳細は過去の予測区間の評価指標の記事を参照してください。

class PredictionIntervalCoverageProbability(PIMetricWithTarget):
    def score(
        self,
        y_true: NumericNDArray,
        y_pred: NumericNDArray,
    ) -> NumericNDArray:
        TargetVariableValidator()(y_true)
        PredictionIntervalValidator()(y_pred)
        return ((y_true >= y_pred[:, 0].reshape(-1)) & (y_true <= y_pred[:, 1].reshape(-1))).mean()  # type: ignore


class NormalizedMeanPredictionIntervalWidth(PIMetricWithoutTarget):
    def __init__(self, R: float) -> None:
        self.R = R

    def score(
        self,
        y_pred: NumericNDArray,
    ) -> NumericNDArray:
        PredictionIntervalValidator()(y_pred)
        return (y_pred[:, 1] - y_pred[:, 0]).reshape(-1).mean() / self.R  # type: ignore

LightGBMを用いたノンパラメトリック分位点回帰による区間予測

　GBDTを用いた分位点回帰は、機械学習によるテーブルデータの区間予測においてベースラインとなるべき標準的なアプローチであると筆者は考えています。例えば90%予測区間を算出したい場合は、5%分位点を予測するモデルと95%分位点を予測するモデルの2つのモデルを学習し、それぞれの予測結果を予測区間とします。特徴量の値に応じた条件付き分位点を予測することになるので、モデルが予測結果に自信を持てるインスタンスでは区間の幅が短くなり、自信を持てないインスタンスでは区間の幅が広くなるような、不均一分散性を考慮した予測区間を得ることが出来ます。

　 LightGBM や XGBoost を用いれば、質的変数や欠損値を含む特徴量でも特段の前処理なく予測モデリング出来ます。さらに、学習時のパラメータを設定するだけでとても簡単に分位点回帰を実装出来るのも大きな長所です。

学習と予測

　今回の実装では、クロスバリデーションを用いて out-of-fold prediction を算出し、これを予測値として評価に使います*1。2つの分位点回帰モデルの学習の実装は下記の通りです:

significance = 0.1
params = {
    "objective": "quantile",
    "learning_rate": 0.01,
    "bagging_fraction": 0.5,
    "feature_fraction": 0.5,
}

y_pred_upper, y_pred_lower = np.empty_like(y), np.empty_like(y)
k_fold = KFold(n_splits=5, shuffle=True, random_state=1)
for train_index, val_index in tqdm(k_fold.split(X)):
    train_X_fold, val_X_fold = X[train_index], X[val_index]
    train_y_fold, val_y_fold = y[train_index], y[val_index]

    train_dataset = lgb.Dataset(train_X_fold, train_y_fold)
    val_dataset = lgb.Dataset(val_X_fold, val_y_fold, reference=train_dataset)
    pred_dataset = lgb.Dataset(val_X_fold, val_y_fold)

    upper_bst = lgb.train(
        {**params, **{"alpha": 1 - significance / 2}},
        train_dataset,
        num_boost_round=10000,
        valid_sets=[val_dataset],
        callbacks=[lgb.early_stopping(100)],
    )

    lower_bst = lgb.train(
        {**params, **{"alpha": significance / 2}},
        train_dataset,
        num_boost_round=10000,
        valid_sets=[val_dataset],
        callbacks=[lgb.early_stopping(100)],
    )

    y_pred_upper[val_index] = upper_bst.predict(pred_dataset.data)
    y_pred_lower[val_index] = lower_bst.predict(pred_dataset.data)

2つのモデルを個別に学習するこのアプローチでは、一部のインスタンスで予測区間の下限が上限を超えてしまう場合があるため、この点を考慮して区間を修正する後処理を行います。

ill_defined_pred = y_pred_upper - y_pred_lower < 0
y_pred_upper[ill_defined_pred] = y_pred_lower[ill_defined_pred]
y_pred_region = np.array([y_pred_lower, y_pred_upper]).T

評価

PICP

PICPを測ってみましょう。90%予測区間を予測出来るよう学習しましたが、実際のPICPは下記の通り83%となりました:

PredictionIntervalCoverageProbability().score(y, y_pred_region)
<<< np.float64(0.8324612403100775)

NMPIW

次に、NMPIWを算出します。正規化に使う定数 の値を指定する必要がありますが、ここでは他の研究でよく用いられている目的変数の最大と最小の差を採用します。実際の の値は約50万ドルとなりました:

y.max() - y.min()
<<< np.float64(4.85002)

NMPIWの値は約25%となり、平均的な区間の幅は目的変数の値の範囲の25%相当であることを示しています:

NormalizedMeanPredictionIntervalWidth(R=y.max() - y.min()).score(y_pred_region)
<<< np.float64(0.24615052721362943)

NMPIWはあくまで平均的な区間の幅に関する指標です。不均一分散を考慮した予測区間では、インスタンスごとに区間の幅は異なります。定数 で除算する前の区間の幅 (10万ドル単位) の分布を可視化してみましょう:

fig = px.histogram(
    x=y_pred_upper - y_pred_lower, title="Prediction Interval Width の分布"
)
fig.update_xaxes(title_text="Prediction Interval Width")

図 (1): 分位点回帰による Prediction Interval Width の分布

一部のインスタンスで幅が大きくなるような、非対称な分布であることが読み取れます。中央値は約10万ドル未満で、平均値よりも小さくなっています。区間の幅が広すぎるとユーザー体験に悪影響があるような場合、(正規化した)区間の幅の四分位点を評価指標としてトラッキングするのも良いでしょう。

Conformalized Quantile Regression

　分位点回帰による区間予測では、元々90%と想定していた予測区間の有意水準に対して、実際のPICPは83%と有意水準を下回りました。本セクションでは、CQRを用いて90%のPICPを達成出来るように分位点回帰による区間を補正してみます。CQRの詳しい解説は過去記事の機械学習による区間予測入門 ④: Conformalized Quantile Regressionを参照してください。

　CQRのアルゴリズム自体はとても簡単です。学習時には分位点回帰による予測区間から下記の通り定数Qを学習します。予測時には分位点回帰による予測区間の上限と下限をQだけ広げてあげるのみです。

E = np.max(np.array([y_pred_region[:, 0] - y, y - y_pred_region[:, 1]]), axis=0)
Q = np.quantile(E, (1 - significance) * (1 + 1 / y.shape[0])).item()
Q
<<< 0.07107777295093937

今回の機械学習タスクではQ = 0.07なので、分位点回帰による予測区間の上限と下限をそれぞれ一律に約7,000ドル広げてあげれば、90%のPICPを達成出来ることになります。下記のセクションでは、これをどのような設計で実装すべきかを考えてみます。

設計

MAPIEとその欠点

　まずは既存のフレームワークで解決出来ないかを考えてみます。scikit-learn-contribプロジェクトのひとつであるとMAPIEというフレームワークを使うのがスタンダードであると思われます。scikit-learn-contribとは、scikit-learnと互換性のある高品質な機械学習フレームワークを収録したGitHub Organizationのことです。imbalanced-learn などもその一員ですが、使ったことがある方も多いのではないでしょうか。

mapie.readthedocs.io

　MAPIEは広く使われているOSSですが、実装上の欠点もあります。CQRのもとになる分位点回帰モデルと、CQR本体の実装が密結合になっている点です。LightGBMを用いたCQRのチュートリアルのページを見てみましょう。CQRの学習と予測が可能なMapieQuantileRegressorに、scikit-learn 互換の学習器をコンストラクタインジェクション出来る作りになっています。分位点回帰にLightGBMを用いる場合は、これにLGBMRegressorを渡すことが出来るのです。この点はクライアント(ユーザー)向けのインターフェースとして良い面もある一方、CQRが持つ model agnostic な長所を活かしきれていない実装になってしまっています。本来のCQRは分位点回帰モデルの予測値さえあれば良く、モデルオブジェクト自体不要なはずです。例えば今回のようにLightGBMのBoosterオブジェクトを学習して得た OOF prediction を元にCQRの学習をしたいようなケースでは、MAPIEを使うことは出来ません。本記事では練習も兼ねて自前の設計を考えてみましょう。

ドメイン駆動で考える機械学習アルゴリズム

ドメイン駆動設計とは

　ソフトウェア工学では、ドメイン駆動設計(DDD)と呼ばれる設計の考え方があります。Wikipediaのドメイン駆動設計のページやドメインモデルのページでは、以下のように定義されています:

ドメイン駆動設計（ドメインくどうせっけい、英語: domain-driven design、DDD）は主要なソフトウェア設計手法の一つであり、ドメインエキスパートの言葉に基づき、ドメインにおけるプロセスやルールをよく表現したドメインモデルを構築し、それに基づいてソフトウェア開発を行うことに主眼を置くものである。

ドメインモデル（英: Domain model）は、システムに関わるさまざまな実体とそれらの関係を説明するシステムの概念モデルである。

　では機械学習システム、とりわけ機械学習アルゴリズムにおけるドメインモデルとは何でしょうか。これは筆者の考えですが、機械学習アルゴリズムが持つ性質や周辺の手法との位置づけを整理したものと言えるのではないでしょうか。これをコードに反映させることがMLシステムの開発で重要なことで、そのためにはアルゴリズムについて理解するだけでなく、サーベイをしっかり読んで周辺の手法との関連や分野全体について知ることが重要だと考えます。上述の意味でのドメイン駆動な設計で機械学習アルゴリズムを実装出来れば、機械学習の研究者が考えるのと同じように手法を追加する (実装を拡張する) ことが容易になるはずです。

CQRのドメインモデルと実装

　では、CQRのドメインモデルを整理してみましょう。CQRとは、Conformal Prediction と呼ばれる予測区間を算出する手法群のひとつです。CPに共通する特徴は、対象とする機械学習タスクに関して model agnostic であることと言えます。対象となる機械学習タスクの予測値を算出し、これを正解データと照らし合わせて予測区間を算出する、というアプローチです。従って、CPは目的変数の予測値と目的変数との関係を学習するスタッキング学習器であると解釈することが出来ます。すなわち、CPは教師あり学習の手法のひとつであると言えるでしょう。この関係性を図示すると以下のようになります:

図 (2): Conformalized Quantile Regression に関するドメインモデルのイメージ。筆者作成。

　では (ソフトウェア設計の観点での) 教師あり学習とは一体何でしょうか。ここからは実装寄りの話になるのですが、scikit-learn の実装を出発点に考えてみると、fitとpredictの両方のメソッドを持ったインターフェースのことであると言えるでしょう。scikit-learn は機械学習実装の標準として優れたインターフェースを提供してくれている一方、いくつか設計上の欠点も抱えていると筆者は考えています。例えば以下の通りです:

• 学習と予測の両方の責務を負っている
• fitは学習済みパラメータを更新するミュータブルな処理であるにも関わらず、自分自身のインスタンスを返す
• データのバリデーションは主にプライベートメソッドの中で行われている

これらの点を克服出来るように、scikit-learn ライクな教師あり学習のインターフェースを整理してみたのが下記の実装です:

class TrainableParams(BaseModel): ...


class BaseSurpervisedLearner[T: TrainableParams](ABC):
    @abstractmethod
    def fit(
        self,
        X: NumericNDArray,
        y: NumericNDArray,
        trainable_params: T | None = None,
    ) -> T:
        raise NotImplementedError


class BasePredictor[T: TrainableParams](ABC):
    @abstractmethod
    def predict(self, X: NumericNDArray, trainable_params: T) -> NumericNDArray:
        raise NotImplementedError


class BaseSupervisedModel(ABC):
    @abstractmethod
    def fit(self, X: NumericNDArray, y: NumericNDArray) -> None:
        raise NotImplementedError

    @abstractmethod
    def predict(self, X: NumericNDArray) -> NumericNDArray:
        raise NotImplementedError


class ValidatableSupervisedModel[T: TrainableParams](BaseSupervisedModel):
    def __init__(
        self,
        learner: BaseSurpervisedLearner[T],
        predictor: BasePredictor[T],
        X_validator: BaseNDArrayValidator,
        y_validator: BaseNDArrayValidator,
    ) -> None:
        self._learner = learner
        self._predictor = predictor
        self._X_validator = X_validator
        self._y_validator = y_validator
        self._trained_params: T | None = None

    @property
    def trained_params(self) -> T | None:
        return self._trained_params

    @property
    def is_trained(self) -> bool:
        return self.trained_params is not None

    def _validiate_is_trained(self) -> None:
        if not self.is_trained:
            raise ValueError("The model is not trained yet.")

    def fit(
        self,
        X: NumericNDArray,
        y: NumericNDArray,
    ) -> None:
        self._X_validator(X)
        self._y_validator(y)
        self._trained_params = self._learner.fit(X, y)

    def predict(self, X: NumericNDArray) -> NumericNDArray:
        self._validiate_is_trained()
        self._X_validator(X)
        return self._predictor.predict(X, self.trained_params) # type: ignore

ValidatableSupervisedModelクラスのコンストラクタに注入する以下の部品を組み替えることで、個別の機械学習アルゴリズムを表現するようなイメージです。

• learner: 学習アルゴリズム
• predictor: 予測アルゴリズム
• X_validator: 特徴量のバリデーター
• y_validator: 目的変数のバリデーター

内部的に学習と予測の責務を分けていること、データのバリデーターをDI出来るようにしたこと、learnerとpredictorが持つメソッドがイミュータブルな処理であることなどが scikit-learn と比較したときの大きな特徴です。

　教師あり学習が抽象であるとすればCPはその具体という位置づけになります。CPは有意水準 alpha をコンストラクタの引数に持ち、スタッキング学習器であるため fit と predict メソッドの引数の名前をわかりやすく呼び替えるようにします:

class BaseConformalPredictionRegressor[T: TrainableParams](ValidatableSupervisedModel[T]):
    def __init__(
        self,
        alpha: float,
        learner: BaseSurpervisedLearner[T],
        predictor: BasePredictor[T],
        X_validator: BaseNDArrayValidator,
        y_validator: BaseNDArrayValidator,
    ) -> None:
        self.alpha = alpha
        super().__init__(learner, predictor, X_validator, y_validator)

    def fit(
        self,
        y_pred: NumericNDArray,
        y_true: NumericNDArray,
    ) -> None:
        return super().fit(X=y_pred, y=y_true)

    def predict(self, y_pred: NumericNDArray) -> NumericNDArray:
        return super().predict(X=y_pred)

　最後に、CQR の実装に移ります。学習と予測のアルゴリズムは前述の通りで、ConformalizedQuantileRegressor クラスがクライアント (ユーザー) 向けの実装となります:

class CqrTrainableParams(TrainableParams):
    Q: float


class CqrLearner(BaseSurpervisedLearner[CqrTrainableParams]):
    def __init__(self, alpha: float) -> None:
        self.alpha = alpha

    def fit(
        self,
        y_pred: NumericNDArray,
        y_true: NumericNDArray,
        trainable_params: CqrTrainableParams | None = None,
    ) -> CqrTrainableParams:
        E = np.max(np.array([y_pred[:, 0] - y_true, y_true - y_pred[:, 1]]), axis=0)
        return CqrTrainableParams(
            Q=np.quantile(E, (1 - self.alpha) * (1 + 1 / y_true.shape[0])).item()
        )


class CqrPreditor(BasePredictor[CqrTrainableParams]):
    def predict(
        self,
        y_pred: NumericNDArray,
        trainable_params: CqrTrainableParams,
    ) -> NumericNDArray:
        return np.array(
            [y_pred[:, 0] - trainable_params.Q, y_pred[:, 1] + trainable_params.Q]
        ).T


class ConformalizedQuantileRegressor(BaseConformalPredictionRegressor):
    def __init__(self, alpha: float) -> None:
        super().__init__(
            alpha,
            CqrLearner(alpha),
            CqrPreditor(),
            PredictionIntervalValidator(),
            TargetVariableValidator(),
        )

　設計の意図が見えやすい実装になっているのではないでしょうか。次のセクションでは実際に ConformalizedQuantileRegressor クラスを用いて California Housing dataset データによる学習と精度評価をしてみます。

CQRの学習と予測

　LightGBMを用いたノンパラメトリック分位点回帰による区間予測のセクションの実装では、クロスバリデーションを用いて out-of-fold prediction を得ました。ここではこの OOF prediction の結果である y_pred_region を用いて CQR のモデルを学習します。評価用のホールドアウトは用意せず、CQRの学習に用いたデータセットで予測し精度評価します。CQRの記事で述べた通り、CQRの学習は low variance な処理なのでリークによる精度の過大評価の影響がほとんどないためです。

cqr = ConformalizedQuantileRegressor(alpha=significance)
cqr.fit(y_pred_region, y)
y_pred_cqr = cqr.predict(y_pred_region)

　分位点回帰では有意水準を90%に設定したにもかかわらず、PICPが83%という結果となりました。CQRによりPICPを90%に近い (well calibrated な) 予測区間にすることが出来ます:

PredictionIntervalCoverageProbability().score(y, y_pred_cqr)
<<< np.float64(0.9)

　CQRは予測区間を定数で補正する手法なので、今回のケースでは区間の幅が広くなっているはずです。下記の通り、分位点回帰のNMPIWは約25%でしたが、CQRにより約28%まで増加しました:

NormalizedMeanPredictionIntervalWidth(R=y.max() - y.min()).score(y_pred_cqr)
<<< np.float64(0.2754608281818479)

　区間の幅の平均的な水準を示すNMPIWだけでなく、データ全体での区間の幅の分布の変化を見てみましょう。CQRは幅を定数で変化させるため、下記グラフのように分布がシフトします:

pi_width = pl.DataFrame(
    {
        "分位点回帰": y_pred_upper - y_pred_lower,
        "CQR": y_pred_cqr[:, 1] - y_pred_cqr[:, 0],
    }
).unpivot(on=["分位点回帰", "CQR"], variable_name="予測器", value_name="pi_width")
fig = px.histogram(
    pi_width,
    x="pi_width",
    color="予測器",
    title="Prediction Interval Width の分布",
    barmode="overlay",
)
fig.update_xaxes(title_text="Prediction Interval Width")

図 (3): 分位点回帰とCQRによる Prediction Interval Width の分布

Calibration set の有無による精度評価への影響

　前セクションでは、リークの影響が小さいと考えられるため Calibration set (CQR用のホールドアウト) を用いずに精度評価しました。本セクションではこのリークの影響による精度の過大評価が本当に生じないのかを検証してみます。ここではCQRの学習をクロスバリデーションを用いて行い、CQRの予測値も OOF prediction とします*2。これによりPICPの低下やNMPIWの増加がみられれば、リークの影響があるため本来は calibration set を用いるべきであることを示しています。

y_pred_cqr_oof = np.empty_like(y_pred_cqr)
k_fold = KFold(n_splits=5, shuffle=True, random_state=1)
cqrs = []
for train_index, val_index in tqdm(k_fold.split(X)):
    train_X_fold, val_X_fold = y_pred_region[train_index], y_pred_region[val_index]
    train_y_fold, val_y_fold = y[train_index], y[val_index]

    cqr = ConformalizedQuantileRegressor(alpha=significance)
    cqr.fit(train_X_fold, train_y_fold)
    y_pred_cqr_oof[val_index] = cqr.predict(val_X_fold)
    cqrs.append(cqr)

foldごとの学習済みパラメータQは、いずれも約0.07と分散が小さくなっています:

[round(c.trained_params.Q, 3) for c in cqrs]
<<< [0.073, 0.07, 0.07, 0.071, 0.072]

PICPとNMPIWを測ってみると、calibration set を用いないケースと精度に差がないことがわかります:

PredictionIntervalCoverageProbability().score(y, y_pred_cqr_oof)
<<< np.float64(0.9002422480620155)

NormalizedMeanPredictionIntervalWidth(R=y.max() - y.min()).score(y_pred_cqr_oof)
<<< np.float64(0.27547253634400953)

今回の California Housing dataset を含めた多くの機械学習タスクでは calibration set 無しでも問題無いと筆者は考えます。ただし、有意水準が厳しいケースや学習データのサンプルサイズが小さいケースではQの分散が大きくなるため過学習が起こり得ることに注意が必要です。

*1:Early stopping に用いているホールドアウトでの予測値をOOF predictionとして評価指標の算出に使っています。この分位点回帰のタスクではモデル選択と評価を同じホールドアウトで行っても精度の過大評価はほとんどみられませんでした。この点はCQRによる精度向上という本記事の主題とあまり関係ないため実験結果を省略します。

*2:本来は calibration set を用いるケースと用いないケースとでサンプルサイズを統制したり、ブートストラッピングにより評価指標の経験分布を求めて指標の平均値の差の検定などをする必要がありますが、ここでは単純な比較で十分であるため行いません。

機械学習による区間予測入門 ⑤: 予測区間の評価指標

　本記事は機械学習による区間予測入門シリーズの5記事目です。他記事へのリンクはこちら: 

機械学習による区間予測入門 ①: はじめに

機械学習による区間予測入門 ②: 分位点回帰

機械学習による区間予測入門 ③: Conformal Prediction

機械学習による区間予測入門 ④: Conformalized Quantile Regression

機械学習による区間予測入門 ⑥: 実装例

目次:

はじめに

　これまでの記事を通じて、様々な機械学習ドメインで分位点回帰モデルを用いた予測区間の算出が機能するケースが多いこと、有意水準を満たす予測区間に補正するアプローチとして Conformalized Quantile Regression がよく用いられていることを紹介してきました。本記事では、機械学習を用いて予測区間を算出する際に、予測区間をどのような指標で評価出来るかを簡単に紹介します。基本的には、「区間の中にどれくらい正解があるか」という精度に関する指標と、「区間は十分に狭いか」を示す広さに関する指標が主です。なお、予測区間の評価指標単体について体系的にまとまったサーベイは (筆者の知る限りでは) 存在せず、ニューラルネットワークを用いた区間予測の手法の提案という文脈で扱われていることがあるのみ、という印象です。これらのペーパーはレファレンスをご参照ください。

評価指標

Prediction Interval Coverage Probability

　PICP は、以下の式 (1) で表される通り、正解ラベルが予測区間の中に入っている確率を表します。予測区間の精度を示す最も基本的な評価指標です。

\begin{align} PICP = \dfrac{1}{N}\sum_{i=1}^N \mathbb{I}(y_i \in [L_i, U_i]). \end{align}

はインジケータ関数、 はそれぞれ予測区間の上限と下限です。

Normalized Mean Prediction Interval Width

　区間の上限と下限の差を、評価データ全体についての平均値に集計したものを Mean Prediction Interval Width と呼びます:

\begin{align} MPIW = \dfrac{1}{N}\sum_{i=1}^N (L_i - U_i). \end{align}

MPIW は平均的な区間の幅を表しますが、目的変数のスケールによって指標のスケールも変わってしまい、異なるドメイン同士での指標の比較が難しくなってしまいます。そこで、目的変数の範囲を示す何らかの定数 で MPIW をスケールした指標を Normalized Mean Prediction Interval Width と呼びます。

\begin{align} NMPIW = \dfrac{1}{N}\sum_{i=1}^N \dfrac{(L_i - U_i)}{R}. \end{align}

としてよく用いられるのは、評価データの目的変数の最大値と最小値の差である です。NMPIW は、区間の広さを表す指標として、PICP と並んでよく用いられてます。R の計算を除けば、正解ラベルが無くとも推論結果のみから測れる指標である点が特徴的です。

　NMPIW ではスケールされた区間の広さの指標を測ることが出来ますが、このスケーリングが不向きなケースもあります。例えば、絵画の資産価値の区間予測問題を考えます。10万円の絵画を±1万円の範囲で予測するのと、1億円の絵画を±1千万円で予測するのとでは、後者の方が区間が広いと考えるべきでしょうか。「±何パーセントの区間なのか」が重要なこのような問題設定の場合には、 の代わりに目的変数 、または回帰モデルによる条件付き期待値 (または分位点回帰モデルによる条件付き中央値) の予測値 をインスタンス毎に除算する下記のような指標が適しているかもしれません:

\begin{align} \dfrac{1}{N}\sum_{i=1}^N \dfrac{(L_i - U_i)}{\hat{y}_i}. \end{align}

ではなく で除算する場合には、正解データが区間からどれくらい外れているかまで含めて測ることになり、純粋な広さの指標と呼べなくなるケースもあるでしょう。他方、 で除算する場合には の予測誤差というノイズが乗ることになります。モデルによる予測区間の予測の外し方と の予測の外し方にある程度連動がある場合や回帰モデルを回しやすい場合は、後者が良いかと思います。なお、式 (4) の指標は筆者が実務において使用したことがありますが、筆者の知る限りではこのような指標を用いている研究は見つけられませんでした。

Mean Absolute Deviation from Prediction Interval

　こちらも筆者が独自に利用したことはあるが他のペーパーでは見たことがない指標になります。④の記事中の non-conformity measure である を、予測を外している企業についてのみ集計した指標になります。

\begin{align} \dfrac{1}{N_{\text{error}}}\sum_{i=1}^{N_{\text{error}}} \max \{L_i - y_i, y_i - U_i \} \end{align}

は予測を外していたインスタンスの数です。予測を外していたインスタンスが、平均的にどれくらい区間から遠くに外していたかを示しています。この指標も目的変数の点予測値や目的変数自身でスケーリング可能です:

\begin{align} \dfrac{1}{N_{\text{error}}}\sum_{i=1}^{N_{\text{error}}} \dfrac{\max \{L_i - y_i, y_i - U_i \}}{\hat{y}_i} \end{align}

また、予測の外し方の外れ値の影響が大きい場合は標本平均ではなく標本中央値や四分位点で測ることも出来ます。

　本記事までの連載を通じて、区間予測の手法と評価指標について触れることが出来ました。最後となる次記事では、実際のデータを使って CQR による予測区間の算出と評価の実装例を紹介します。

レファレンス

• Cordier, T., Blot, V., Lacombe, L., Morzadec, T., Capitaine, A., & Brunel, N. (2023, August). Flexible and Systematic Uncertainty Estimation with Conformal Prediction via the MAPIE library. In Conformal and Probabilistic Prediction with Applications (pp. 549-581). PMLR.
• Khosravi, A., Nahavandi, S., Creighton, D., & Atiya, A. F. (2010). Lower upper bound estimation method for construction of neural network-based prediction intervals. IEEE transactions on neural networks, 22(3), 337-346.
• Khosravi, A., Nahavandi, S., & Creighton, D. (2010). Construction of optimal prediction intervals for load forecasting problems. IEEE Transactions on Power Systems, 25(3), 1496-1503.
• Khosravi, A., Nahavandi, S., Creighton, D., & Atiya, A. F. (2011). Comprehensive review of neural network-based prediction intervals and new advances. IEEE Transactions on neural networks, 22(9), 1341-1356.

機械学習による区間予測入門 ④: Conformalized Quantile Regression

　本記事は機械学習による区間予測入門シリーズの4記事目です。他記事へのリンクはこちら: 

機械学習による区間予測入門 ①: はじめに

機械学習による区間予測入門 ②: 分位点回帰

機械学習による区間予測入門 ③: Conformal Prediction

機械学習による区間予測入門 ⑤: 予測区間の評価指標

機械学習による区間予測入門 ⑥: 実装例

目次:

はじめに

　分位点回帰のみを用いる予測区間の算出方法だと不均一分散性を考慮した柔軟な予測区間を作成できる一方、事前に決めた有意水準を満たせるかどうかが非自明であったり、実証的にも有意水準を満たせない場合があることを②の記事で確認しました。他方、 ③の記事の通り、Conformal Prediction はモデル非依存な予測区間を算出でき、予測区間の中に正解ラベルが含まれる確率について望ましい性質を持っていますが、区間の幅が定数になってしまう制約がありました。本記事で紹介する Romano et al. (2019) による Conformalized Quantile Regression は、両者の長所を組み合わせたアプローチとして予測区間の計算で広く応用されています。

	分位点回帰	CP	CQR
予測区間の幅	◯	×	◯
予測区間の性質	△	◯	◯

　本記事は目次の通り、アプローチの直感的理解や具体的なアルゴリズムを先に説明し、その後でアルゴリズムの解釈や実務上のポイント、手法の限界について述べます。また、最後のオプショナルなセクションでは well-calibrated property と呼ばれる有限標本性質を証明していますが、議論の大枠は③の記事に基づいています。

CQR の直感的理解

　CQRのエッセンスは大きく2点です:

1. 分位点回帰モデルを学習して予測区間の上限と下限を予測する
2. 1のモデルの予測結果をどれくらい補正したら良いかを、ホールドアウトしておいたデータ (calibration set) で計算する

　1では分位点を予測出来るモデルであれば好きな学習アルゴリズムを用いてモデルを学習して良いです。例えば LightGBM を用いて上限と下限の2つのモデルを学習したり、ニューラルなモデルであればモデルを2つ学習するやり方の他に区間そのものを学習する方法もあります*1。ここが大きな利点であり、欠損値や質的変数を含む特徴量から成るテーブルデータのMLタスクの場合は LightGBM で分位点回帰が可能で、自然言語や画像認識などのMLタスクはそのドメインに合ったニューラルネットワークの損失関数を pinball loss に設定することで分位点回帰を行えるので、非常に広範な機械学習ドメインで予測区間を計算することが出来る手法なのです。

　2のステップでは、分位点回帰による予測結果が予め定めた予測区間の有意水準を満たしていない場合、CQRによる区間を広げる補正が可能です。逆に分位点回帰による予測区間が不必要に広い場合は、正解データが区間の中に入る確率を元の有意水準に保ったまま区間を狭めるよう補正出来ます。2のステップは1で予測した予測区間と正解データにのみ依存し、1の学習アルゴリズムには依存せずに区間の補正が可能です。

図 (1): CQRで補正したい予測区間の例。

アルゴリズム

ここでは前セクションの1,2のステップの具体的な手順である、アルゴリズムの詳細をまず紹介します。

1. 予測区間の有意水準 を定める
2. データを train set と calibration set に分割する
3. Train set を用いて、有意水準が となるように分位点回帰モデルを学習する
4. Calibration set を用いて、3のモデルによる各インスタンスの予測結果 を得る
5. Calibration set を用いて、式 (1) で表される各インスタンスの non-conformity measure を計算する

\begin{align} E_i = \max \{\hat{q}_{\text{lower}}(x_i) - y_i, y_i - \hat{q}_{\text{upper}}(x_i)\} \end{align}

1. Calibration set における について、有意水準 の empirical quantile である定数 を計算する
2. 未知のインスタンス に対する予測区間は次の通り算出する:

疑似コードは以下のようなイメージです:

# 有意水準を定める (90%予測区間なら alpha = 0.1)
alpha = 1 - 0.9

# 分位点回帰モデルの学習は好きなやり方で良い
# ここでは5%分位点と95%分位点を予測する2つのモデルを学習する
model_lower = train_lower_quantile_reg(X_train, y_train, significance=alpha / 2)
model_upper = train_upper_quantile_reg(X_train, y_train, significance=1 - alpha/2)

# calibration set における予測値を得る
y_cal_lower = model_lower.predict(X_cal)
y_cal_upper = model_lower.predict(X_cal)

# Non-conformity measure である E を計算し、E をもとに定数 Q を算出する
E = np.max(np.array([y_cal_lower - y_cal, y_cal - y_cal_upper]).T, axis=1)
Q = np.quantile(E, ((1 - alpha) * (1 + 1 / X_cal.shape[0])))

# 未知データへの予測にはモデルの予測値に ± Q の定数で補正する 
y_test_lower = model_lower.predict(X_test) - Q
y_test_upper = model_upper.predict(X_test) + Q

アルゴリズムの解釈

　下記の証明のセクションで示す通り、未知のインスタンス におけるこのアルゴリズムで導かれた 予測区間である は、 "well calibrated property" と呼ばれる以下の性質を満たします:

\begin{align} 1 - \alpha \leq P(y_{\text{test}} \in C(x_{\text{test}})) \leq 1 - \alpha + \dfrac{1}{n_{\text{cal}}+1} \end{align}

　ここでは上記のアルゴリズムのセクションの5〜7のステップがどのような意味を持つのか、何を計算してどのように式 (2) を達成しているのかを考えてみます。

Non-conformity measure

　式 (1) における calibration set 内のあるインスタンス に関して、モデルの予測値である を固定値として、正解ラベル が変化したときに non-conformity measure がどのように変化するかを図示したのが図 (2) です。

図 (2): あるインスタンス の正解ラベル と non-conformity measure の関係。

は予測区間が正解とどれだけ不一致しているかを示す尺度で、予測区間の中に正解ラベルがある場合は値が負になり、予測区間の外にある場合は正、区間の端にちょうど正解が一致した場合はゼロとなります。例えばアイスクリームの価格を100円〜120円と予測していたときに、実際の正解が130円だった場合、 円だけ区間の端から予測を外している、という意味です。正解ラベルや予測区間そのものとスケール (単位) が一致した尺度であることがポイントで、これにより後述の予測区間の補正が可能になります。

予測区間の補正

　前述のような を calibration set の全てのインスタンスで計算してみたときに、 はどのような分布をしているでしょうか。例えばアイスクリームの価格予測の例で 予測区間を求めたいケースを考えます。 の約 分位点を集計した結果、 が得られたとします。このとき、calibration set 全体のうちの約 が となっていて、残りの約 が となります。そのため分位点回帰モデルの予測結果のみで を達成でき、区間を補正する必要はありません。

　次に、上述のアイスクリームの価格予測の 予測区間の例で、 円が得られたケースを考えます。 このとき、 なので、 となっているインスタンスは全体の に満たないはずです。分位点回帰による予測区間が狭すぎて、 となるインスタンスが多い状態です。すなわち分位点回帰モデルの予測結果のみだと となり、区間を広くするよう補正してあげる必要があります。 の約 分位点では区間の左右どちらかの端から正解データまで 円分予測を外してしまっているので、未知のインスタンスの予測時に一律に区間の上限と下限の両方を 円分広くしてあげれば、 を達成出来るはずです。

　逆に負の値である 円が得られた場合はどうでしょうか。 なので、 となっているインスタンスは全体の より多い状況です。分位点回帰モデルの予測結果のみだと区間が広すぎて となり、もっと区間を狭くしても を維持出来ます。 の約 分位点では 円だけ区間が左右どちらかに余分に広い状況なので、未知のインスタンスの予測時に区間の上限と下限の両方を 円分狭くしても に出来るはずです。

　下の図 (3) は区間が広すぎるケースと狭すぎるケースにおける、分位点回帰による予測区間と の分布を示しています。どういう状況で がゼロから乖離するのか、イメージが湧いてくると思います。

図 (3): 分位点回帰による予測区間と non-conformity measure の分布の例。左は なので CQR により区間は狭くなる補正がかかり、右は なので CQR により広くなる。

実務上の注意点

ホールドアウト　

　CQR のアルゴリズムに厳密に従うならば、分位点回帰モデルのモデル選択に用いる validation set とは別に、 の集計に用いる calibration set という新たなホールドアウトが必要となります。ですが筆者の経験上、モデル学習時の validation set と CQR の calibration set は同じデータを使っても問題無いケースが多いと思います。 は に関するシンプルな piecewise linear な変換を行っているだけで、 はその分位点の集計値にすぎず、非常に low variance な変換となります。そのため同じホールドアウトを二重に使ってもリークによる精度の過大評価の影響はほとんど無いケースが大半だと思います。もちろん、元の分位点回帰モデルをクロスバリデーションにより out-of-fold prediction を出力する形で学習し、CQR を OOF prediction に対するスタッキング学習器であると解釈すれば、ホールドアウトを増やすかどうかで悩む必要は無いでしょう。他方、巨大なニューラルネットワークで学習に時間を要するなど、気軽にクロスバリデーションを行えないケースはあります。その場合は、上記の理由で calibration set を別にホールドアウトしない、すなわち validation set を使い回すという選択肢もあり得ると思います。

有意水準

　分位点回帰のみではなくCQRを用いることの追加的な利点として、有意水準を固定することで問題設定をシンプルにしやすいという点もあげられます。例えば、予測区間算出におけるユーザー体験として望ましい有意水準 が既知であるとします。このとき、CQR を用いずに分位点回帰のみで区間予測のモデリングをすることを考えてみます。未知データにおける予測区間で を満たせなかった場合、特徴量エンジニアリングといった分位点回帰モデルのモデル選択に加えて、 を別の値である にチューニングすることで元の を達成出来ないかも検討しなければなりません。他方、CQR を用いる場合は、 は元の値に固定した上で分位点回帰モデルのパフォーマンスを如何に向上させるかに集中することが出来ます。なぜなら、 は CQR による区間の補正で達成出来るからです。もちろん、望ましい有意水準 自体が探索の対象であることの方が多いでしょうが、その場合でも望ましい有意水準の探索の問題と有意水準が固定となる分位点回帰モデルのパフォーマンス向上の問題とを分けて考えることが出来るようになります。

CQR の限界

予測区間の周辺確率と条件付き確率

　 ③の記事で解説した Conformal Prediction と同様に、式 (2) はあくまで について周辺化した確率についての主張で、 がどのような値でも成り立つ訳ではありません。すなわち、以下の式 (3) のような条件付き確率に関する命題が成り立つ訳ではありません。

\begin{align} 1 - \alpha \leq P(y_{\text{test}} \in C(x_{\text{test}}) | x_{\text{test}} = x ) \leq 1 - \alpha + \dfrac{1}{n_{\text{cal}}+1}. \end{align}

CQRは分位点回帰の予測結果を定数 で補正しているだけなので、通常は図 (4) の真ん中の図のように精度が高いところもあれば低いところもあってデータ全体で集計すると式 (2) が成り立っている状態となります。

図 (4): 予測区間の命題が周辺確率について成り立っていることを示す例 (中央)。Angelopoulos & Bates (2022) より転載。

分位点回帰の重要性

　CQR を使えば、元の分位点回帰モデルがどんなものであっても、どれほど精度が悪くても 式 (2) を満たすような補正が可能です。一方、CQR は前セクションの通り式 (2) の周辺確率を満たすために定数 で予測区間を補正するアプローチなので、そもそも元の条件付き分位点回帰の予測精度がある程度高いことが重要になります。裏を返せば、CQR 無しの分位点回帰のみでも十分に精度が高い状態でないと扱いづらいという欠点を抱えた手法でもあります。分位点回帰の精度が低い場合は、図 (5) の右側のように による補正の影響が大きくなり、全体的に区間の幅が広いような予測結果になってしまいます。

図 (5): 精度の良い分位点回帰 (左) と精度の悪い分位点回帰 (右) による CQR の違い。

　このように、予測区間の算出においては、「区間の広さ」といった 以外の観点も重要になります。次回の記事では、機械学習における予測区間の評価指標について紹介します。

Well-calibrated property (オプショナル)

　本セクションでは、CQR による予測区間の有意水準が有限標本においても満たされることを示す、式 (2) を証明します。

　最も重要な仮定として、 Conformal Prediction の場合と同様に、確率変数列 の exchangeability を仮定します。Exchangeability は IID をより一般化したような性質で、ほとんどの機械学習タスクにおいて満たされる性質と考えて差し支えないです。詳細が気になる方は ③の記事を参照してください。

　 を、以下の通り未知データにおける non-conformity measure とします:

\begin{align} E_{\text{test}} : = \max \{\hat{q}_{\text{lower}}(x_{\text{test}}) - y_{\text{test}}, y_{\text{test}} - \hat{q}_{\text{upper}}(x_{\text{test}})\} \end{align}

CQR による予測区間の定義から、

\begin{align} y_{\text{test}} \in C(x_{\text{test}}) \iff E_{\text{test}} \leq Q_{1–\alpha}(E, \mathbb{D}_{\text{cal}}). \end{align}

従って、

\begin{align} P(y_{\text{test}} \in C(x_{\text{test}})) = P(E_{\text{test}} \leq Q_{1–\alpha}(E, \mathbb{D}_{\text{cal}})). \end{align}

さらに、 および を算出する学習モデルが、③の記事で定義した permutation invariant であるとします*2。 が exchangeable であるため、 for も exchangeable となります。従って、③の記事の命題 1 と同様の証明により、

\begin{align} 1 - \alpha \leq P(E_{\text{test}} \leq Q_{1–\alpha}(E, \mathbb{D}_{\text{cal}})) = P(y_{\text{test}} \in C(x_{\text{test}})) \leq 1 - \alpha + \dfrac{1}{n_{\text{cal}}+1}. \end{align}

レファレンス

• Angelopoulos, A. N., & Bates, S. (2021). A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511.
• Kuchibhotla, A. K. (2020). Exchangeability, conformal prediction, and rank tests. arXiv preprint arXiv:2005.06095.
• Romano, Y., Patterson, E., & Candes, E. (2019). Conformalized quantile regression. Advances in neural information processing systems, 32.

*1:詳細は②の記事の脚注 1 を参照してください。

*2:元論文では触れられていませんが、 の exchangeability を導くには、学習モデルの permutation invariance の仮定が必要になります。③の記事の脚注 6 で述べた通り、勾配ブースティングやニューラルネットワーク等のバギングやミニバッチを用いる学習アルゴリズムは厳密には permutation invariant ではありませんが、 exchangeability preserving であるとみなせるため、permutation invariance の仮定のみでも一般性を失うことはありません。実務上のインプリケーションとしては、学習データの順序が重要になるような例外的な学習アルゴリズムでなければ、予測区間の補正に non-conformity measure の分位点を用いるアプローチが機能する、ということです。詳細は該当の脚注を参照してください。

機械学習による区間予測入門 ③: Conformal Prediction

　本記事は機械学習による区間予測入門シリーズの3記事目です。他記事へのリンクはこちら: 

機械学習による区間予測入門 ①: はじめに

機械学習による区間予測入門 ②: 分位点回帰

機械学習による区間予測入門 ④: Conformalized Quantile Regression

機械学習による区間予測入門 ⑤: 予測区間の評価指標

機械学習による区間予測入門 ⑥: 実装例

目次:

はじめに

　Conformal Prediction とは、データの分布や学習アルゴリズムに依存せずに、有限標本においても統計的に望ましい性質を満たした予測区間を算出出来るアルゴリズムのことです。このアルゴリズムの解析的な性質の研究や拡張・応用を含めた広義の分野を指して Conformal Prediction と呼ぶこともあり、初学者向けのチュートリアルやサーベイ論文も多く出版されています (Angelopoulos & Bates, 2021; Fontana et al., 2023; Lei et al., 2018; Shafer & Vovk, 2008)。他方、このアルゴリズムの統計的性質の厳密な導出を一部スキップしてしまっているものが多く、本記事はそのギャップを埋める self contained な記述と直感的な分かりやすさを両立することを目指しました。本記事は前記事の分位点回帰とは独立した内容になっていますが、分位点回帰とCPを理解することが後続の Conformalized Quantile Regression の記事の内容の理解の助けになると思います。尚、本記事をスキップしても、次記事の内容を大枠で理解する上では支障ありません。

　本記事の構成は以下の通りです。最初に、CP の具体的なアルゴリズムを含む概要について紹介します。次に、CP の性質を導出する上で最も重要な仮定である、確率変数の exchageability の定義とその性質を述べます。次に、前セクションの性質を用いて、CP が持つ有限標本性質を導出します。最後に、CP の限界について触れます。なお、定義や証明のパートは主に Kuchibhotla (2020) と Romano, et al (2019) を参考に再構成しています。

Conformal Prediction の概要

　本セクションではCPのアルゴリズムとその性質を紹介します。厳密なことより前に具体的なアルゴリズムを知ることで、やろうとしていることの直感やその性質に対するイメージを持ってもらうことがゴールです。

アルゴリズム

　CPは解くべきMLタスクが分類・回帰であるかを問わず、またどのような学習アルゴリズムであっても予測区間を算出出来る手法です。分析者は、 "non-conformity measure" と呼ばれる、モデルの予測値が正解ラベルの実測値とどれくらい離れているかを示すヒューリスティックなスコアを定義する必要があります。CPはヒューリスティックな non-conformity measure を変換して厳密な予測区間に変換してくれる、図 (1) のようなイメージです。このプロセスは学習データとは別にホールドアウトした calibration set と呼ばれるデータセットで行います。学習データで回帰や分類モデルを学習し、calibration set ではモデルの予測誤差に基づいて区間を計算します。具体的には次のような手順のアルゴリズムです:

1. データを train set と calibration set に分割する
2. train set でモデルを学習する
3. 予測区間の有意水準 を定める
4. calibration set で、モデルの予測値と正解ラベルの実測値に基づいて non-conformity measure を算出する
5. calibration set で、non-conformity measure の 分位点 を算出する
6. non-conformity measure が 以下となるような 目的変数の範囲を予測区間とする

図 (1): ヒューリスティックな non-conformity measure を厳密な予測区間に変換する CP のアルゴリズム。Anastasios et al. (2022) より転載。

回帰の例

　例えば回帰の場合、non-conformity measure にはL1誤差 (Mean Absolute Deviation) を用いることが多いです。95%予測区間を算出する場合の実装例は以下のようになります:

model.fit(X_train, y_trian)

alpha = 1 - 0.95
scores = np.abs(y_cal - model(X_cal))
significance = np.ceil((n_cal + 1) * (1 - alpha)) / n_cal
q_hat = np.quantile(scores, significance, method="higher")

# 未知データの予測区間を算出する (L1誤差 q_hat から逆算する)
prediction_interval_upper = model(X_new) + q_hat
prediction_interval_lower = model(X_new) - q_hat

CPの性質の直感的理解

Well calibrated な性質

　上記のプロセスで計算された予測区間 は、「回帰モデルの split conformal prediction」のセクションで示す通り、予測区間が有意水準 とほぼ等しくなるような "well calibrated property" と呼ばれる下記の有限標本性質を満たします:

\begin{align} 1 - \alpha \leq P(Y_{\text{new}} \in \hat{C_{n_{\text{cal}}}^{\alpha}}) \leq 1 - \alpha + \dfrac{1}{n_\text{cal}+1}. \end{align}

予測モデルの精度の重要性

　モデルの予測性能が悪いと区間はどうなるでしょうか。scores の推定値がどれも大きい値のものばかりになり、結果 q_hat の値が大きくなり、予測区間の幅も大きくなってしまいます。そのため、CPはモデルの精度に依らずに有意水準を満たす予測区間を計算出来ますが、十分に狭い予測区間を算出するには元のモデルの精度が高いことが重要です。

Non-conformity measure の単調変換

　non-conformity measure をL1誤差ではなくL2誤差 (y_cal - model(X_cal)) ** 2 で測った場合はどうでしょうか。「回帰モデルの split conformal prediction」のセクションで述べる通り、CP で算出する予測区間は non-conformity measure の単調変換に対して不変となります。直感的には、CPは「学習に使っていない未知データの誤差を有意水準 % の分位点に収めるには、予測値をどれくらい補正したら良いか」を計算しているアルゴリズムであるため、誤差の分位点 (順位) のみが重要になるのです。

定数となる予測区間の幅

　CPによる予測区間は、モデルによる予測値に対して定数 2 * q_hat となる予測区間を算出します。そのため、図 (2) のように目的変数が不均一分散性を持つような分布に従っている場合でも、特徴量の値を条件づけたときの区間を予測の不確実性に応じて広くしたり狭くしたりすることが出来ません。

図 (2): 不均一分散性を持つ目的変数の分布に対する Conformal Prediction。Fontana, et al. (2023) より転載。

　ここまでで具体的なアルゴリズムを交えてCPの概要に触れました。ここからはCPとは何か、どのような仮定の下でどのような性質を導けるかを追っていきます。次セクションでは、CPにおいて最も中心的な役割を果たす確率変数の Exchangeability の定義とその性質について紹介します。その次のセクションでは、Exchageability のセクションで導出した諸定理を用いてCPとその性質を導出します。

Exchangeable な確率変数とその性質

Exchangeability の定義

　 とします。実数の確率変数 () について、 から までのいかなる順列 でこの確率変数を並べ替えても、並び替える前と後との同時確率分布が等しいとき、この確率変数は exchangeable であると言います。すなわち、任意の について下記が成り立つことが exchangeability の定義です:

\begin{align} P(Y_1 \leq a_1, ..., Y_n \leq a_n) = P(Y_{\pi(1)} \leq a_1, ..., Y_{\pi(n)} \leq a_n). \end{align}

　例として、下記の同時正規分布に従う確率変数 は exchangeable です。

\begin{align} \begin{bmatrix} X_1 \\ X_2 \end{bmatrix} \sim N \begin{pmatrix} \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \end{pmatrix}. \end{align}

例えば だと一方が高い値なら他方も高い値な傾向にあるので互いに独立ではありませんが、 と を入れ替えても同じ同時分布に従っているので exchangeable であると言えます。

　Exchangeability は IID を一般化したような性質で、 exchangeable な確率変数列は identical な分布に従っている必要がありますが、必ずしも独立である必要はありません*1。IID の分布に従っていることは exchageability の十分条件ではあるが、必要条件ではありません。

Exchangeable な確率変数と順位

　ここでは、exchangeable な確率変数とその順位に関する定理とその系について紹介します。

定理 1

実数の確率変数 が exchangeable であるとき、これらの順位は の分布に依らず の全ての順列の数だけ一様に生じ得ます: \begin{align} (\text{rank}(Y_i; \{Y_1, ..., Y_n\}): i \in [n]) \sim \text{Unif}(\{\pi: [n] \rightarrow [n]\}). \end{align}

ここで、 は 以下の値を取る要素の数を表す関数です。また、 は の全ての順列に関する一様分布を表しており、各順列は の確率で生じ得ることを意味しています。

証明:

ここでは証明の簡略化のために が almost surely distinct (値に重複が無い、同順位な値が無い確率が1) であると仮定します*2。この仮定により、事象 は全ての順列 について互いに素となります。また、いずれかひとつの組み合わせが起こり得るので、 \begin{align} \sum_{\pi: [n] \rightarrow [n]} P(\{Y_1 \leq ... \leq Y_n \}) = 1. \end{align} ここで、 とすると、exchangeability の仮定より*3、

\begin{align} P(Y_1 \leq ... \leq Y_n) = & P( (Y_1, ..., Y_n) \in A_n) \\ = & P( (Y_{\pi(1)}, ..., Y_{\pi(n)})\in A_{\pi(n)}) \\ = & P(Y_{\pi(1)} \leq ... \leq Y_{\pi(n)}). \end{align}

従って、式 (5) より となります。順位の定義から、 ] が の特定の順列となる事象は と同値であるため、定理 1 が証明されます。

　定理 1 から次の系を導くことが出来ます:

系 1

を flooring (小数を切り捨てて整数を得る) 関数とすると、 \begin{align} P\left(\text{rank}(Y_i; \{Y_1, ..., Y_n\}) \leq t\right) = \dfrac{\lfloor t \rfloor}{n}. \end{align} さらに はP値として有効な値として次の性質を満たします: \begin{align} P(p \leq \alpha) \leq \alpha \quad \forall \alpha \in [0, 1]. \end{align}

証明:

\begin{align} P(\text{rank}(Y_i; \{Y_1, ..., Y_n\}) \leq t) &= P(\text{rank}(Y_i; \{Y_1, ..., Y_n\}) \leq \lfloor t \rfloor) \\ &= \sum_{j=1}^{\lfloor t \rfloor} P(\text{rank}(Y_i; \{Y_1, ..., Y_n\}) = j) \\ &= \sum_{j=1}^{\lfloor t \rfloor} \dfrac{(n-1)!}{n!} \end{align}

ここで とすると、 となります。

Exchangeable な確率変数の変換

定理 2

Exchangeable な確率変数 について、確率変数 への変換 を考えます。全ての順列 について、 \begin{align} \pi_m G(x) = G(\pi_n x) \quad \forall x \in \mathbb{X}^n \end{align} を満たすような順列 が存在することと、 が exchangeability preserving である ( が exchangeable であるときに が exchangeable である) ことは同値となります。

証明:

同時確率分布が等しいことを示す記号を とします。全ての について式 (14) を満たす が存在するとすると、 のexchangeability より、

\begin{align} G(X) & \stackrel{d}{=} G(\pi_n X), \quad \forall \pi_n \\ & \stackrel{d}{=} G(\pi^*_n X) \\ & \stackrel{d}{=} \pi_m G(X), \quad \forall \pi_m. \end{align}

従って が exchangeability preserving となります。逆の場合は本記事の後続の議論では扱わないため証明を割愛します。

Conformal Prediction

　本セクションでは、Conformal Prediction による予測区間とその性質を導出します。確率変数列の exchangeability を仮定するのみで導くことが出来るのがこの手法の特徴的な点です。第一に、最も単純なケースとして、 exchangeable な確率変数列 について の予測区間を求める、単変数である確率変数の Full conformal prediction について扱います。第二に、同様に単変数である確率変数の Split conformal prediction と呼ばれるアプローチを扱います。目的変数だけでなく特徴量のデータが得られるケースへの拡張に必要なステップです。最後に、第二のアプローチの拡張として、 のような特徴量と目的変数のペアである確率変数列が得られる場合に回帰モデルを学習した上で予測区間を計算する split conformal prediction について扱います。

単変数の Full Conformal Prediction

　Conformal prediction のエッセンスを理解するため、まずは最も単純なケースとして、 exchangeable な確率変数 を得られる場合の予測区間について考えます。 を ceiling 関数とし、 とすると、] であるので、系 1 の式 (9) より下記が成り立ちます:

\begin{align} P\left(\text{rank}(Y_{n+1}; \{Y_1, ..., Y_{n+1}\}) \leq \lceil (n+1)(1-\alpha)\rceil \right) = \dfrac{\lceil (n+1)(1-\alpha) \rceil}{n+1}. \end{align}

さらに、下記の通り予測区間 を定義します:

\begin{align} \hat{C_n^{\alpha}} := \{y \in \mathbb{R}: \text{rank}(y; \{Y_1, ..., Y_n, y\}) \leq \lceil (n+1)(1-\alpha)\rceil \} \end{align}

ここで、 を に関する有意水準 の empirical quantile とします。分位点の定義より、予測区間は以下のように empirical quantile を用いて表すことが出来ます:

\begin{align} \hat{C_n^{\alpha}} = \{y \in \mathbb{R}: y \leq \hat{Q}_{n+1}(1-\alpha) \}. \end{align}

Empirical qunatile が 個の確率変数列について計算されているのは、 に関する順位に基づいているためです。この empirical quantile は、 を除いた のみから計算することも出来るので、そのような変形を試みます。まず、 を、 の中の 番目に小さい値とします。すると、任意の について、以下が成り立ちます*4:

\begin{align} Y_{n+1} \leq Y_{(k,n)} \iff Y_{n+1} \leq Y_{(k,n+1)}. \end{align}

ここで、分位点の定義より、 かつ、 であるため、

\begin{align} Y_{n+1} \leq \hat{Q}_{n+1} (1-\alpha ) \iff Y_{n+1} \leq \hat{Q}_{n} ( (1-\alpha ) (1 + 1/n ) ). \end{align}

従って、

\begin{align} \hat{C_n^{\alpha}} = \{y \in \mathbb{R}: y \leq \hat{Q}_{n} ( (1-\alpha ) (1 + 1/n ) ) \}. \end{align}

予測区間を のみのempirical qunatile を用いて表すことが出来ました。ここで、 であり、 であるため、確率変数列 における の予測区間 について、"well calibrated property" と呼ばれる下記の命題が成り立ちます:

命題 1

確率変数列 が exchangeable であるとすると、 \begin{align} 1 - \alpha \leq P(Y_{n+1} \in \hat{C_n^{\alpha}}) \leq 1 - \alpha + \dfrac{1}{n+1}, \quad \forall \alpha \in [0, 1], n \geq 1. \end{align}

　このアプローチを特徴量と目的変数の両方の確率変数列 が得られる場合に拡張していくための次のステップとして、次のセクションでは引き続き単変数の確率変数列が得られるケースを考えます。

単変数の Split Conformal Prediction

　 の exchangeable な確率変数列があるとします。ここでは、 の平均値を予測値とする単純な学習モデルを構築出来る場合に、どのように予測区間を計算出来るかを考えてみます。 を train set と calibration set の2つに分割します。train set で平均値 を計算し、calibration set で non-conformity measure を計算します。 下記の式 (25) のように、calibration set 内の non-conformity measure の順位、あるいは empirical quantile に基づいて 予測区間 を定義します*5:

\begin{align} \hat{C_{n_{\text{cal}}}^{\alpha}} : & = \{y \in \mathbb{R}: \text{rank}(r; \{r_1, ..., r_{n_{\text{cal}}}, r\}) \leq \lceil (n_{\text{cal}}+1)(1-\alpha)\rceil \} \\ & = \{y \in \mathbb{R}: r \leq \hat{Q}_{n_{\text{cal}}} ( (1-\alpha ) (1 + 1/n_{\text{cal}} ) ) \}. \end{align}

は定理 2 より明らかに exchangeability preserving となり、命題 1 と同様の証明で下記の命題が成り立ちます:

命題 2

確率変数列 が exchangeable であるとすると、 \begin{align} 1 - \alpha \leq P(Y_{n+1} \in \hat{C_{n_{\text{cal}}}^{\alpha}}) \leq 1 - \alpha + \dfrac{1}{n_\text{cal}+1}, \quad \forall \alpha \in [0, 1], n \geq 1. \end{align}

　本セクションでは、目的変数のみを用いてその平均値を予測する、単純な回帰モデルを学習出来る場合の conformal prediction について扱いました。次セクションでは一般の回帰モデルの場合について扱います。 Exchangeable な特徴量から予測値への非線形な写像である学習モデルが、exchangeability preserving となることがカギとなります。

回帰モデルの Split Conformal Prediction

　Exchangeable な特徴量と目的変数のペア が得られるとします。再び train set と calibration set にデータを分割することを考えてみます。学習データ を用いて学習済みモデル を得ることが出来ます。 Calibration set においてモデルの予測値と正解ラベルの実測値とのL1 誤差で測った を non-conformity measure とします。このとき、 は exchangeable となるかを考えてみましょう。

まず、関数 が次に定義されるような permutation invariant な性質を満たすとします:

\begin{align} g(z; z_1, ..., z_n, z_{n+1}) = g(z; z_{\pi(1)}, z_{\pi(n)}, z_{\pi(n+1)}) \end{align}

このとき、式 (14) を満たす が存在するため、 は exchangeability preserving です。 が permutation invariant な学習モデルである場合*6、non-conformity measure も exchangeability preserving となります。従って、上セクションと同様に train set でモデルを学習して calibration set で non-conformity measure を計算することで、下記の性質を満たす予測区間を形成出来ます:

命題 3

確率変数列 が exchangeable であるとすると、 \begin{align} 1 - \alpha \leq P(y_{n+1} \in \hat{C_{n_{\text{cal}}}^{\alpha}}) \leq 1 - \alpha + \dfrac{1}{n_\text{cal}+1}, \quad \forall \alpha \in [0, 1], n \geq 1. \end{align}

　ここでは non-conformity measure として L1 誤差を用いましたが、 L2 誤差など単調変換となる non-conformity measure では、値の順序が変わらないため予測区間は等しくなります。元のモデルの精度が良いか悪いかにかかわらず式 (29) の性質は満たされますが、精度が悪いと区間の計算に必要な non-conformity measure の分位点の値が大きくなるため、その分区間の幅が広くなってしまいます。

Conformal Prediction の限界

予測区間の周辺確率と条件付き確率

　式 (29) はあくまで calibration set 全体において予測区間が正解を含む周辺確率についての命題であり、calibration set 内の個々のインスタンスについて特徴量を条件づけた条件付き確率が成り立っていることを主張している訳では無いことに注意が必要です。すなわち、下記が成り立つ訳ではありません:

\begin{align} 1 - \alpha \leq P(y_{n+1} \in \hat{C_{n_{\text{cal}}}^{\alpha}} | x_{n+1} ) \leq 1 - \alpha + \dfrac{1}{n_\text{cal}+1}, \quad \forall \alpha \in [0, 1], n \geq 1. \end{align}

通常は図 (2) のように精度が高いセグメントもあれば低いセグメントもあり、calibration set 全体では式 (29) を満たしている状態になります。

予測区間の幅

　「Conformal Prediction の概要」のセクションで述べた通り、CPは区間の幅が定数となってしまう扱いづらさがあります。この欠点を克服したのが、Conformalized Quantile Regression と呼ばれる手法で、分位点回帰とCPの利点を組み合わせたものになります。次の記事では実務上使いやすいこの手法について紹介します。

レファレンス

• Angelopoulos, A. N., & Bates, S. (2021). A gentle introduction to conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.07511.
• Fontana, M., Zeni, G., & Vantini, S. (2023). Conformal prediction: a unified review of theory and new challenges. Bernoulli, 29(1), 1-23.
• Kuchibhotla, A. K. (2020). Exchangeability, conformal prediction, and rank tests. arXiv preprint arXiv:2005.06095.
• Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., & Wasserman, L. (2018). Distribution-free predictive inference for regression. Journal of the American Statistical Association, 113(523), 1094-1111.
• Romano, Y., Patterson, E., & Candes, E. (2019). Conformalized quantile regression. Advances in neural information processing systems, 32.
• Shafer, G., & Vovk, V. (2008). A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3).
• Vovk, V., Gammerman, A., & Shafer, G. (2005). Algorithmic learning in a random world (Vol. 29). New York: Springer.

*1:証明はKuchibhotla (2020) を参照してください。

*2:重複がある (同順位があり得る) 場合でも、一様分布に従う と任意の を考え、 とすると、 exchangeable な確率変数とIIDな確率変数の和は exchangeable であることを示すことができ、almost surely distinct な連続確率変数列 が得られるので、元の証明でも一般性を失いません。 は実務上は1e-08など限りなく小さい値に設定することで、同順位のときにはいずれかをほぼ等確率で選ぶように設定出来る、という意味のパラメータです。

*3: が exchangeable であることは、任意のBorel measurable set について、 が成り立つことと同値でもあります(詳細はKutibhotla (2020)を参照してください)。例えば とすると、 は 閉集合である と の共通集合あるため、Borel measurable setです。同様に はBorel measurable set であるため、exchangeability の定義により式 (6) から (7) への等号が成り立ちます。

*4: であるとき、 となり、必要条件を満たします。十分条件は であるため成り立ちます。

*5:定理 1 の場合と同様に の almost sure distinctness を一般性を失うことなく仮定出来ます。

*6:勾配ブースティングやニューラルネットワークなどの確率的な性質を持つモデルは、厳密には学習データの順序に対して permutation invariant ではありません。例えば、バギングやミニバッチで抽出されるデータの順序が入れ替わると出力されるモデルが変わり得ます。しかし、例えば学習データを並び替えたモデルによる予測値が、元のモデルの予測値に対してIIDのホワイトノイズを加えただけの結果になると仮定しても実用上問題無い場合、すなわち となる場合は、non-conformity measure が exchangeable となるため、命題 3 が同様に成り立ちます。実務上は、学習データの順序が重要なアルゴリズムでない限り、学習モデルは permutation invariant であると考えて差し支えないはずです。

機械学習による区間予測入門 ②: 分位点回帰

　本記事は機械学習による区間予測入門シリーズの2記事目です。他記事へのリンクはこちら: 

機械学習による区間予測入門 ①: はじめに

機械学習による区間予測入門 ③: Conformal Prediction

機械学習による区間予測入門 ④: Conformalized Quantile Regression

機械学習による区間予測入門 ⑤: 予測区間の評価指標

機械学習による区間予測入門 ⑥: 実装例

目次:

確率変数の分位点とその推定方法

分位点とは

　確率変数 の 分位点とは、直感的に言うと確率 %以上で実現するような の値の最小値のことです。例えば ] の一様分布に従う確率変数の5%分位点は0.05、といった具合です。50%分位点のことを特に中央値と呼びます。25%と75%分位点のことを四分位点とも呼び、データの分布 (の実現値) の性格を表す集計値としてよく見ている方も多いでしょう。厳密には、確率変数 の分布関数を とすると、 の 分位点 は次式の通り定義されます: 

\begin{align}
q_Y(\alpha) = \inf \{y: F_Y(y) \geq \alpha \}
\end{align}

の逆関数が存在するときはシンプルに となるのですが、 が離散値をとる場合など分布関数の逆関数が存在しない場合も含めたより一般的な定義では式 (1) のようになります。下の図 (1) を参照するとinfimumを用いた定義のイメージが湧きやすいかもしれません。

図(1): 確率変数 の分布関数とその分位点の例。WikipediaのQuantile function のページより転載。

次のセクションでは確率変数の分位点がどのように推定出来るかを解説します。

分位点の推定方法

　確率変数の分位点は、 と定義された "pinball loss" と呼ばれる損失関数の期待値の最小化により得られることが知られています ( はインジケータ関数): 

\begin{align}
q_Y(\alpha) = \arg \min_u E[\rho_\alpha(Y - u)]
\end{align}

　まずは式 (2) の右辺の損失関数の意味を直感的に理解しましょう。損失関数は次のようにも表現できます: 

\begin{align}
\rho_\alpha(Y - u) = \begin{cases}
  \alpha (Y - u) & \text{if } Y - u \geq 0 \\    
  -(1-\alpha) (Y - u), & \text{otherwise.}
\end{cases}
\end{align}

、 すなわち中央値の場合は平均絶対誤差となる L1 loss の期待値を最小化することと同値です。下の図 (2) は のケースの損失関数をグラフとして描いたものです。 のように が小さい値の場合、分位点の推定値が大きくなることに大きなペナルティが課され、分位点の推定値が小さくなることには小さなペナルティが課される結果、中央値よりも小さい推定値が得られます。 が大きい値の場合はその逆となります。

図 (2): alpha = 0.1 の pinball loss。筆者作成。

　Pinball lossの直感はざっくりと掴めたでしょうか。次に、式 (2) を証明します。ここでは計算の簡略化のために、 が ] の一様分布に従う場合のみに限定しましょう。分布の仮定に依らず、pinball lossの期待値は次のように表現できます。

\begin{align}
E[\rho_\alpha(Y - u)] = (\alpha - 1) \int_{-\infty}^u (y - u)dF_Y(y) + \alpha \int_u^\infty (y - u)dF_Y(y)
\end{align}

が ] の一様分布に従う場合、

\begin{align}
\min_{u \in [0,1]} E[\rho_\alpha(Y - u)] &= (\alpha - 1) \int_0^u (y - u)dy + \alpha \int_u^1 (y - u)dy \\
&= \dfrac{1}{2}\left\{u^2(1-\alpha) + \alpha(1-u)^2\right\} 
\end{align}

式 (6) の最小化の1階条件を解くと、

\begin{align}
\dfrac{d}{du}\left(E[\rho_\alpha(Y - u)]\right) = 0 \iff u = \alpha
\end{align}

] の一様分布では であるため、  となり、pinball lossの期待値最小化により得られる推定量は 分位点と等しくなります。

線形分位点回帰 

線形分位点回帰 (LQR) とは

　分位点回帰 (quantile regression) とは、特徴量から分位点を予測するモデル、すなわち、特徴量 で条件づけた目的変数 の条件付き 分位点を推定する回帰モデルのことです。パラメータに関して線形なモデルの場合を特に線形分位点回帰 (linear quantile regression) と呼びます。線形モデルのOLS推定により得られる回帰モデルが の条件付き期待値を表す Conditional Expectation Function であるのに対して、分位点回帰では Conditional Quantile Function が得られます。

図 (3): 線形分位点回帰の例。scikit-learnのQuantile regressionのページより転載。

モデルパラメータの推定方法

　線形分位点回帰モデルも上セクションの場合と同様に、目的変数と線形モデルとの L1 誤差をベースとした pinball loss の最小化によりパラメータを推定出来ます (Koenker & Bassett, 1978)。

\begin{align}
\beta_{\alpha} = \arg \min_{\beta \in \mathbb{R}^K} E[\rho_\alpha(y - X\beta)]
\end{align}

推定アルゴリズム

　線形分位点回帰における pinball loss の最小化問題は、線形計画問題の標準形に変形することが出来るため、シンプレックス法と呼ばれる線形計画アルゴリズムで解くことが出来ます。目的変数と予測値との差分が正となる example と負となる example とに分け、またパラメータも正と負の部分とに分けることでLPの形に上手く変形していきます。 まずは、

\begin{align}
\mu &= |y - X\beta|\mathbb{I}(y - X\beta \gt 0) \\
\nu &= |y - X\beta|\left\{ 1-\mathbb{I}(y - X\beta \gt 0)\right\}
\end{align}

とすると、式 (8) の sample counterpart は以下のように変形できます: 

\begin{align}
& \min_{\beta} \imath^T \mu + \imath^T \nu \\
\text{s.t.} \quad & y - X\beta = \mu + \nu \text{, } \{\mu, \nu\} \in \mathbb{R}^N_{+}
\end{align}

ここで は column of ones とします。さらに , とし、以下の通り定義すると、

\begin{align}
\tilde{X} &= [X^T, -X^T, I, -I] \\
\theta &= [\beta_{+}^T, \beta_{-}^T, \mu^T, \nu^T]^T \\
z &= [\mathbf{0}^T,\mathbf{0}^T,\imath^T,\imath^T]^T
\end{align}

pinball loss の最小化問題は以下のように表現出来ます: 

\begin{align}
& \min_{\theta} z^T\theta \\
\text{s.t.} \quad & \tilde{X}\theta = y, \quad \theta \geq 0.
\end{align}

上記の定式化が線形計画問題の標準形となります。scikit-learnのQuantileRegressorの実装を見ると同様にLPで解いていることが分かります。

LQRを用いた予測区間の限界

　では線形分位点回帰モデルを使って予測区間を作ることは出来るでしょうか。例えば5%分位点と95%分位点を予測する2つの線形モデルを学習してその間を90%予測区間とすることは出来ます。しかし、このやり方だと model specification が線形に限られてしまうだけでなく、予測区間の幅もパラメータ推定値に関して線形な値である となってしまいます。「特徴量がこういう値のときには精度が悪いので区間の幅は広くなり、こういう値では狭くなる」といった表現力に乏しい状況です。この点を克服するためには、非線形な条件付き分位点回帰モデルが必要になります。次セクションではこのような選択肢のうちの一つである、GBDTを用いたノンパラメトリック分位点回帰について紹介します。

GBDTによるノンパラメトリック分位点回帰

GBDTによる分位点回帰

　欠損値や質的変数を含む特徴量から成るタビュラーデータなど、GBDTによるモデリングが適したMLタスクである場合、GBDTを用いて分位点回帰モデルを実装するのが現実的な選択肢となります。XGBoost や LightGBM など任意の損失関数をインジェクト可能なGBDTのフレームワークでは、パッケージ内に pinball loss が実装されておりとても手軽に分位点回帰モデルを学習出来るためです。例えばLightGBMでは、次のように train 関数にパラメータを渡すだけです: 

lgb.train(params={"objective": "quantile", "alpha": 0.05, ...}, ...)

この実装で式 (3) の通りの損失関数でGBDTのモデルを学習することが出来ます。

GBDTを用いた予測区間の限界

　線形モデルの場合と同様に、区間の下限と上限を予測する2つのGBDTモデルを学習することで、予測区間を算出することが出来ます*1。筆者の経験上、ベースラインのモデリングとしてはこのやり方で実務上十分に機能するケースが多いという印象です。ですが、XGBoostやLightGBMによる分位点回帰では、予測結果の性質について良く知られていません。他の学習アルゴリズムでも、漸近的な性質が知られているのみで有限標本性質については分かっていないものが多いです。

　未知データにおいても予測区間の性質 を満たせるようにするための簡便なアプローチとしては、2つの分位点回帰モデルのパラメータ をホールドアウトセットでチューニングするというやり方もあるでしょう。しかし、このやり方だと が解釈しづらくなるという点以外にも重要な欠点があります。分位点回帰モデルの予測精度向上という問題と予測区間の精度向上という問題とを同時に扱わなければならない、すなわち、モデルパラメータのチューニングにより両方のゴールを追わなければならない、という問題設定の複雑さを生んでしまうことです。これらの欠点を克服するためのアプローチとして、Conformalized Quantile Regression と呼ばれる手法があります。次稿では、 CQR の考え方の元となっている Conformal Prediction について紹介します。

レファレンス

• Kabir, H. D., Khosravi, A., Hosen, M. A., & Nahavandi, S. (2018). Neural network-based uncertainty quantification: A survey of methodologies and applications. IEEE access, 6, 36218-36234.

• Koenker, R., & Bassett Jr, G. (1978). Regression quantiles. Econometrica: journal of the Econometric Society, 33-50.

• Koenker, R. (2005). Quantile regression. Cambridge Univ Pr.

• Quantile function. (2024, Sep 29). In Wikipedia. https://en.wikipedia.org/wiki/Quantile_function

*1:NLPやCVなど、ニューラルネットワークが適したMLタスクの場合でも、まずは損失関数をpinball lossとして予測区間の下限と上限を予測する2つのモデルを実装するのが、ベースラインとしては適していると思います。ベースラインの次の一手として、NNの性質を活かして予測区間を作るには Khosravi et al., (2018) のサーベイを出発点とすると良いです。