異なるアルゴリズム間のピッチ推定結果の時間ズレを補正して長さを揃えるPython関数

先日の記事内で引用した論文に、時間ズレ補正のアルゴリズムが提案されていたので、それを Python で使えるようにしてみたかった、ということ。

[2602.01727] Voting-based Pitch Estimation with Temporal and Frequential Alignment and Correlation Aware Selection

例えばこう書く。使いやすさを考慮して、ターゲットのピッチ系列長をリファレンスのピッチ系列長に揃える形とした。

import numpy as np

def align_and_pad(
    ref_f0, target_f0, max_shift: int = 5, rpa_threshold_cents: float = 50.0
):
    """Align and pad a target pitch sequence to match a reference in time.

    This function searches for the optimal frame shift that minimizes the
    logarithmic error between valid (voiced) segments of the two pitch
    sequences, compensating for algorithm-specific padding or windowing
    differences. The target is then shifted and padded with NaNs to match
    the exact length of the reference.

    Args:
        ref_f0 (np.ndarray): Reference pitch array in Hz. Unvoiced frames
            should be represented as np.nan or values <= 0.
        target_f0 (np.ndarray): Target pitch array to be aligned.
        max_shift (int, optional): Maximum number of frames to search for
            the optimal temporal discrepancy. Defaults to 5.
        rpa_threshold_cents (float, optional): Tolerance threshold in
            cents for calculating Raw Pitch Accuracy (RPA). Defaults to
            50.0.

    Returns:
        tuple: A tuple containing:
            - aligned_f0 (np.ndarray): The corrected target pitch array,
              having the exact same length as `ref_f0`.
            - best_shift (int): The optimal shift amount in frames applied
              to the target. Positive means delayed.

    """
    ref_len = len(ref_f0)

    # 1. Calculate the cross-correlation in valid (non-NaN/0) intervals
    # to find the discrepancies.
    # (For simplicity, we find the shift that minimizes the absolute
    # error in frames where both values <200b><200b>exist.)
    best_shift = 0
    max_rpa = -1.0  # RPA is maximized

    for shift in range(-max_shift, max_shift + 1):
        correct_frames = 0
        valid_frames = 0

        for i, tgt_f0 in enumerate(target_f0):
            ref_idx = i + shift

            if 0 <= ref_idx < ref_len:
                v_ref = ref_f0[ref_idx]
                v_tgt = tgt_f0

                # Check if both are valid voiced frames
                if (
                    v_ref > 50
                    and v_tgt > 50
                    and not np.isnan(v_ref)
                    and not np.isnan(v_tgt)
                ):
                    valid_frames += 1

                    # Eq (12): Calculate the deviation in the cent
                    # region
                    diff_cents = np.abs(1200 * np.log2(v_ref / v_tgt))

                    # Eq (11): If it is within the threshold, it will be
                    # counted as "correct answer (1)".
                    if diff_cents < rpa_threshold_cents:
                        correct_frames += 1

        # RPA calculation (only if valid frames exist)
        if valid_frames > 0:
            current_rpa = correct_frames / valid_frames

            if current_rpa > max_rpa:
                max_rpa = current_rpa
                best_shift = shift

    # Use the found `best_shift` to pad and even out the length.
    aligned_f0 = np.full(ref_len, np.nan)
    for i, tgt_f0 in enumerate(target_f0):
        ref_idx = i + best_shift
        if 0 <= ref_idx < ref_len:
            aligned_f0[ref_idx] = tgt_f0

    # aligned_f0 will have exactly the same length as ref_f0 (the
    # missing part is padded with NaN)
    return aligned_f0, best_shift

Pythonで使える基本周波数(ピッチ)抽出ツール

並べてみるとこんな感じ:

抽出方法(CREPE と FCNF0++ を除く):

import librosa
import parselmouth
import pyreaper
import pysptk
import pyworld as pw

filename = "./hoge.wav"

sr = 48000 # Hz
frame_length = 2048 # point
frame_period = 5  # msec
hop_length = int(fs * (frame_period / 1000)) # point
fmin, fmax = librosa.note_to_hz("C3"), librosa.note_to_hz("C5")

y, _ = librosa.load(filename, sr=sr)

f0_rapt = pysptk.sptk.rapt(y.astype(np.float32), sr, hop_length, fmin, fmax)
f0_swipe = pysptk.sptk.swipe(y.astype(np.float32), sr, hop_length, fmin, fmax)
f0_pyin, _, _ = librosa.pyin(
    y=y, fmin=fmin, fmax=fmax, sr=sr, frame_length=frame_length, hop_length=hop_length
)
f0_dio, _ = pw.dio(y.astype(np.float64), sr, frame_period=frame_period)
f0_harvest, _ = pw.harvest(y.astype(np.float64), sr, frame_period=frame_period)
f0_praat = parselmouth.Sound(y, sr).to_pitch().selected_array["frequency"]

x = (np.iinfo(np.int16).max * y).astype(np.int16)
_, _, _, f0_reaper, _ = pyreaper.reaper(
    x, sr, fmin, fmax, frame_period=frame_period / 1000
)

CREPE と FCNF0++ は MacOS 26.4.1, Python 3.14 ではうまくインストールできなかった.Tensorflow が Python 3.14 非対応という点は大きい.

基本周波数系列の推定結果を凸最適化によりポストフィルタする話

はじめに

CA社の記事があった.

AI Lab、音声・音響信号処理分野のトップカンファレンス「ICASSP 2026」にて2本の論文採択 | 株式会社サイバーエージェント

同記事では下記論文が採択された旨が書いてあった.

[2602.01727] Voting-based Pitch Estimation with Temporal and Frequential Alignment and Correlation Aware Selection

基本周波数系列を推定する手法は数多く提案されている.それら手法の推定結果の median を取るヒューリスティックな後処理(多数決,voting)によって,推定誤りを補完し,頑健な推定が可能となることが知られている. 本論文では,voting により推定精度が改善する根拠および voting の妥当性を, median や mode に着目して理論的に示しつつ,さらなる精度改善の手法を提案している.

わたしは複数の推定結果の voting という後処理、いわばポストフィルタに興味を惹かれた.特に median による統合処理には凸最適化が使えそうと分かり,ちょっと手元で推定を試してみたのであった. 本記事はその記録である.

定式化

まず多数の推定器の結果をアンサンブルするために,基本周波数の推定に関する voting では median を,また 有声/無声判定の推定では mode を取っている.その根拠・妥当性は論文に書かれているので詳細は割愛する(サボり).一応,簡易メモを残しておくと,

  • median(基本周波数):推定の誤差分散を具体的に導出し,相関の低い推定器を取り込むことで誤差分散が  1/n オーダー( n は手法数)で減少することを示している.相関の高い推定器を取り込む場合(=どの推定器も同じ間違いをする)は,誤差分散は単独の手法( n = 1)と変わらない.

  • mode(有声/無声判定):二値分類問題として扱われるが,mode による統合はコンドルセの陪審定理(Condorcet's Jury Theorem)によって正当化される.多数決が正解する確率を具体的に示し,推定器の個数  n が増えるにつれて,その確率が 1 に収束することを述べている.

論文ではさらに推定結果を向上させるために「時間シフト」と「周波数シフト」を適用して誤差分布を中心化し,補正済みのデータに対して median を取ることで,誤差分散の縮小を狙っている.

本記事は,研究の発展を見すえて,これらシフトを適用する前段階の medianを取る処理 (基本周波数が対象)を,凸最適化に置き換える試みである.

各手法の推定結果に対する median の計算は,L1ノルム(差の絶対値の和)の最小化に等しい(後述の「補題」参照).凸最適化問題として定式化すると,次式となる.


\begin{align*}
\hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left[ \sum_{t=1}^{T} \sum_{i=1}^{n} | y_{i,t} - \theta_t | + \lambda \sum_{t=2}^{T} | \theta_t - \theta_{t-1} | \right]
\end{align*}

 y_{it} は手法  i による時刻  t の観測値(推定値), \theta_{t} は推定対象となる時刻  t の基本周波数, n は手法数, T はフレーム数である. \lambda は重みパラメタである. 系列全体をベクトルを使って  \boldsymbol{\theta} = (\theta_{1}, \theta_{2}, \ldots, \theta_{T} )^{\top} と書いている.

第1項は観測誤差をL1ノルム(絶対誤差の和)で測っており,推定値を観測データの median に近づける効果がある.データの忠実性(fidelity)を評価しているとも見なせる.また第2項は正則化項であり,フレーム間におけるピッチの変化量を測りつつ,総和を取ることで,フレーム全体の総変動(total variation)を評価している.この項の導入により,系列の不要な振動を抑えて平坦にしつつ,急激な変化も許容する効果を持つ.

ちなみにベクトル記法に移れば,次式となる.


\begin{align*}
\hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}}  \left( \sum_{i=1}^{n} \| \mathbf{y}_i - \boldsymbol{\theta} \|_1 + \lambda \| \mathbf{D} \boldsymbol{\theta} \|_1 \right)
\end{align*}

 \| \cdot \|_1 はL1ノルムである.また  \mathbf{D} は差分行列である:


\begin{align*}
 \mathbf{D} \boldsymbol{\theta} = (\theta_2 - \theta_1, \ldots, \theta_T - \theta_{T-1})^{\top}
\end{align*}

補題

median を取る処理と,L1ノルムの最小化が等しいことを示しておく.

推定したい値を  \theta,観測値を  x_1, \ldots, x_n とする. 目的関数を「L1ノルム」(絶対誤差の和)と置く.


\begin{align*}
J (\theta) = \sum_{i=1}^{n} |x_i - \theta|
\end{align*}

この  J(\theta) を最小にする  \theta を求めたい.そこで  \theta で微分してみる.


\begin{align*}
\frac{d }{d\theta} J (\theta) = \sum_{i=1}^{n} \frac{d }{d\theta}  |x_i - \theta| = \sum_{i=1}^{n} - \textrm{sign}(x_i - \theta)
\end{align*}

ここで 絶対値  |x| の微分が符号関数  \textrm{sign} (x) になることを使った( x \gt 0 なら +1 x \lt 0 なら -1). 上記の導関数を  =0 と置いてみる.


\begin{align*}
\sum_{i=1}^{n} \textrm{sign}(x_i - \theta) = 0
\end{align*}

つまり最適な  \theta は,「 x_{i} \gt \theta となるデータの個数」(プラスの符号)と 「 x_{i} \lt \theta となるデータの個数」(マイナスの符号)が釣り合うような位置に相当する.これは median の定義そのものであり,まさに分布の「中央」に位置しているというわけである.

なぜ L1 ノルムが適切なのか?

基本周波数(ピッチ)の推定において,L2(平均値, mean)ではなくL1(中央値, median)を選ぶ理由は何だろうか. つまるところ、「外れ値への頑健性」にあるといえる.

ピッチ推定では倍ピッチ・半ピッチと呼ばれる推定誤りがしばしば発生する.平均値を計算すると,これら推定誤りの影響を受けやすい. 例えば,各推定器の出力が  [100, 102, 200] (Hz) だとすると,平均値は 134 (Hz) になり,倍ピッチ誤り(= 外れ値)200 (Hz) の影響を強く受けている.これは L2 の明確な弱点である.それに対して中央値は 102 (Hz) になり,外れ値の影響を受けづらい.

実験

実験その1

実験用にダミーデータを生成し,凸最適化が期待通りに動作しているかを確認する.

実験条件

ダミーデータ(真のピッチ)の仕様:

  • 総フレーム数 200,実区間は  [0, 4 \pi ]

  • 第50フレームまでは  100 + 10 \sin(t)

  • 第51フレーム〜第100フレームは  150 (定数)→ジャンプを表現

  • 第101フレーム〜第150フレームは  150 + 20 (t - t_{100}) / 10 → t_{100} はほぼ  2\pi

  • 第151フレーム〜第200フレームは  120 (定数)

観測データ(つまり各推定器の推定値)は真のピッチにノイズを加える.具体的にはラプラスノイズ(loc=0.0, scale=2.0)を加えつつ, オクターブエラーを表現するため 10 %の確率で +1 オクターブ(2倍)もしくは -1 オクターブ(0.5倍)のミスをするようにノイズ(外れ値)を加えた. 以下にデータ生成の実装を抜粋したものを示す.

# 1. 真のピッチ軌跡
t = np.linspace(0, 4 * np.pi, n_frames)
true_pitch = np.zeros(n_frames)
true_pitch[:50] = 100 + 10 * np.sin(t[:50])  # Vibrato
true_pitch[50:100] = 150  # Jump
true_pitch[100:150] = 150 + 20 * (t[100:150] - t[100]) / 10  # Glissando
true_pitch[150:] = 120  # Jump down

# 2. 各手法による推定値
observations = []
for i in range(n_methods):
    # 基本ノイズ (ラプラス分布:裾が重い)
    noise = np.random.laplace(loc=0.0, scale=2.0, size=n_frames)

    # オクターブエラー (確率的に大幅にズレる外れ値を追加)
    # 10%の確率で、+1オクターブ(2倍) or -1オクターブ(0.5倍) のミスをする
    outliers = np.random.choice([1.0, 2.0, 0.5], size=n_frames, p=[0.9, 0.05, 0.05])

    obs = true_pitch * outliers + noise
    observations.append(obs)

観測系列の数 n_methods は 5 とした.

凸最適化の実装には,cvxpy を用いた.実装部分を抜粋する:

import cvxpy as cp

def solve_l1_filtering(observations, lam=5.0):
    n_methods, n_frames = observations.shape

    # 最適化変数 (推定したい真のピッチ)
    x = cp.Variable(n_frames)

    # 1. データ忠実項
    # L1ノルムを使うことで、観測データの中央値(median)に近づこうとする
    # sum_{i} |x - y_i|
    data_fidelity = 0
    for i in range(n_methods):
        data_fidelity += cp.norm(x - observations[i], 1)

    # 2. 正則化項
    # total variation (TV) = 隣接フレームの差分のL1ノルム
    # これにより、急激な変化を許容しつつ、不要な振動を抑える
    regularization = cp.norm(cp.diff(x), 1)

    # 目的関数
    objective = cp.Minimize(data_fidelity + lam * regularization)

    # 制約なし最適化問題として定義
    problem = cp.Problem(objective)

    # ソルバーで解く (ECOSやOSQPなどが自動選択される)
    problem.solve()

    return x.value

lam は正則化の重みである.最適化対象となる目的関数に応じたデータ形式で用意し,あとは Minimize を呼び出せばよいので,とても直感的に使うことができる.

実験結果

図1にピッチ推定の結果を示す.正則化の重みは 5.0 としている.真のピッチ(ground truth, 黒点線)に対して,ノイズありデータ系列はかなり激しく「暴れている」様子が見て取れる.青線が凸最適化を行わず単純に median だけで推定した結果を示し,赤線が今回の凸最適化で推定した結果を示している.median のみでもノイズの影響をそれほど受けずに推定できているが,一部でパルス的な観測値がそのまま残った形になっている.隣接フレームの影響を考慮していないから,といえるだろう.それに対して凸最適化の結果,正則化項の効果により,パルス的な推定値も消え,真のピッチに近い推定値が得られた.また前半部分のビブラート的な波打つ変動にも追従できている.

図1:ピッチ推定の結果(正則化の重み  \lambda = 5.0

図2は正則化項の重みを 100.0 にしたときの結果を示している.スパース性が一層強く出ており,よりノート(音符)らしい推定結果,つまり一部の急激なジャンプ(第50フレーム,第150フレーム)を除いて平坦な系列が得られている.

図2:ピッチ推定の結果(正則化の重み  \lambda = 100.0

実験その2

実際の音声データを対象にして,推定実験を行う.

実験条件

音声データは The LJ Speech Dataset の "LJ001-0001.wav" とした.

ピッチ推定には librosa から使える pYIN のみを使い,異なるピッチ推定器を複数個使わない. ただしピッチの抽出条件を変えて,擬似的な複数の推定器を用意した. 具体的には分析窓の長さを 2048をデフォルトとして,他を1024,4096 としている. 窓を短くすれば(1024)時間分解能が上がるが周波数分解能が下がるので,低音部分の抽出が不安定になる. また長くすれば(4096)安定するが細かい変化には追従しづらい.

実験結果

図3にピッチ推定の結果を示す.正則化の重みは 1.0 とした.スペクトログラムに見える縞模様のうち「一番下」のものをトレースできていれば良い.0.3 〜 0.4 秒付近にかけて,pYIN (win=4096) が推定誤りを起こしていることが分かる.また 1.8 〜 1.9 秒付近では pYIN (win=1096) のみが正しいと思われるピッチ系列を推定しているが,他の推定器では追随できておらず,多数決の結果として平坦な系列が得られている.その他のフレームを見ても,概ね各推定器の結果を統合して滑らかにしたような軌跡が得られている.

図3:ピッチ推定の結果(正則化の重み  \lambda = 1.0

正則化の重みを 10.0 に変えて推定したときの結果を図4に示す.細かな上下動が抑制され,滑らかなピッチ軌跡が得られている.

図4:ピッチ推定の結果(正則化の重み  \lambda = 10.0

最後に,正則化の重みを 50.0 に変えて推定したときの結果を図5に示す.大きな上下動を起こす箇所はわずかであり,ほぼ平坦な軌跡が得られている.これはデメリットではなく,採譜の際にはむしろ有効に働くことが期待される.

図5:ピッチ推定の結果(正則化の重み  \lambda = 50.0

おわりに

凸最適化による基本周波数(ピッチ)推定の後処理はうまく動いたといえる.今回は pYIN の結果のみを統合したが,論文で試されていた他の推定器の結果を統合した実験は実施しておきたいところである(アライメントの問題は適切に対処する前提で).

ピッチ推定手法をさらに発展させるために,以下の方向性が考えられる.

  1. 各推定器の絶対誤差に対して,推定器の信頼度を反映する重みを導入する方法である.重みの決め方には色々なものが考えられる.有声/無声判定にも重みは導入できる(重み付き多数決).

  2. 現状の枠組みは,真値の周りに,単一のラプラス分布に従う誤差が発生しているモデルとして解釈できる.ラプラス分布の尤度を最大化することで,median が得られたのであった. これを拡張するため,ラプラス混合分布を導入する.真値の分布以外に,倍ピッチ・半ピッチ誤りに対応する分布を置く.これにより,推定値が真値から生成される過程をうまくモデル化できることが期待される(生成モデル).

  3. 現状の枠組みは,フレームごとに独立して計算し,後処理的にスムージングをしていた.これを時系列モデルや状態空間モデルに組み込む方向性も考えられる(例:粒子フィルタ,隠れマルコフモデル).論文で議論されていた時間方向・周波数方向のアライメントの問題を,確率モデルによってうまく解決できる可能性がある.

  4. 論文のセクション3.2では「Greedy Algorithm による相関を考慮した手法選択」を提案していた.これは事前の静的な選択であるが,「曲のAメロでは Harvest が良いが、Bメロの高音では pYIN が良い」といった時間変化する信頼度には対応できない.そこで Group Lasso による動的な重み付けを導入する方向性が考えられる.

  5. ピッチの軌跡を,「ゆっくり変化する曲のメロディ」と「ビブラートなどの周期的な揺れ」,そして「突発的なエラー(スパースなノイズ)」に分解して最適化する方向性が考えられる.Robust PCA と類似の凸最適化になる.

  6. 現状の実装は cvxpy のソルバに依存している(今回は ECOS を採用).当然,フレーム数や手法数が増えると最適化に多くの計算時間を要する.そこで ADMM が使えるように問題を再定式化し,高速化を図る方向性が考えられる.実際にそれは可能であって,手元の ADMM 実装を用いた数十秒の歌声データに対する簡易的な実験では,最適化の計算時間を削減しつつ(3割〜4割程度の削減; ハイパパラメタ依存),ほぼ同等のピッチ推定が得られることを確認している.

『統計的学習理論』の第1章のメモ:2値判別問題におけるベイズ規則とベイズ誤差について

p.9 の定義 1.1 を引用する。

損失関数  l を定めたとき、任意の可測関数  h : \mathcal{X} \rightarrow \mathcal{Y} のもとでの予測損失の下限


\begin{align*}
\mathop{\rm inf} \limits_{h: \text{可測関数}} R(h)
\end{align*}

を、損失関数  l のもとでのベイズ誤差(Bayes error)という。下限を達成する仮説が存在するとき、その仮説をベイズ規則(Bayes rule)という。

学習の目標は、学習データからベイズ規則を達成する仮説  h_0 を求める(推定する)ことである。

2値判別問題におけるベイズ規則は p.10 に記述が見つかる。すなわち、「入力  x が与えられたとき、最も出現する確率が大きなラベルを予測ラベルとする仮説が最適」である。そこに書かれている内容で実質的な証明になっているが、もう少し整理しておきたい。

データの確率分布を  D とする。すなわち評価データは  (X, Y) \sim D のように分布するが、 D Y に関して周辺化した分布を  D_{\mathcal{X}} と書く。さらに  X を固定したときの  Y の分布を  D_{Y \mid X} と書く。 また 0-1 損失関数を  l_{\text{err}} と書き、そのときの予測損失(予測判別関数)を  R_{err} とする。 任意の仮説  h: \mathcal{X} \rightarrow \mathcal{Y} について、


\begin{align*}
R_{\text{err}}(h) &= \mathbb{E}_{(X, Y) \sim D} [ l_{\text{err}} (h(X), Y)] \\
&= \mathbb{E}_{X\sim D_{\mathcal{X}}}   [ \mathbb{E}_{Y\sim D_{Y \mid X}} [ l_{\text{err}} (h(X), Y) \mid X] ] \\
&=  \mathbb{E}_{X\sim D_{\mathcal{X}}}   [ \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h(X) \neq Y ]  ] ]
\end{align*}

となる。さらに期待値  \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h(X) \neq Y ]  ] について、


\begin{align*}
 \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h(X) \neq Y ]  ] &= \mathrm{Pr}(Y =1 \mid X)  \mathbf{1}  [ h(X) =0 ]  +  \mathrm{Pr}(Y =0 \mid X)  \mathbf{1}  [ h(X) =1 ]
\end{align*}

と展開できる。ここで


\begin{align*}
\alpha_{X} = \mathrm{min} \{ \mathrm{Pr}(Y =1 \mid X), \mathrm{Pr}(Y =0 \mid X)\}
\end{align*}

と置いてみる。 h(X) =0 のとき期待値は   \mathrm{Pr}(Y =1 \mid X) \geq \alpha_{X} であり、また  h(X) =1 のときも同様に   \mathrm{Pr}(Y =0 \mid X) \geq \alpha_{X} であるから、まとめると以下を得る。


\begin{align*}
 \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h(X) \neq Y ]  ] &\geq  \alpha_{X}
\end{align*}

ゆえに


\begin{align*}
R_{\text{err}}(h) \geq \mathbb{E}_{X\sim D_{\mathcal{X}}}   [   \alpha_{X} ]
\end{align*}

であり、下限(下界)が得られた。

次に下限を達成する仮説が実際に存在することを示したい。式を簡略化するため、 \mathrm{Pr}(Y =1 \mid X) = \beta_{X} とおく。 X=x が与えられたとき、 \beta_{X} \geq 1/2 ならば  1 を出力し、 \beta_{X} \lt 1/2 ならば  0 を出力する仮説  h_{0} を考えてみる。 前者の条件は  \beta_{X} \geq 1- \beta_{X} =  \mathrm{Pr}(Y =0 \mid X = x) と等価であり、後者は  \beta_{X} \lt 1- \beta_{X} と等価である。 両者の条件付き確率が等確率ならば、どちらも 1/2である。 予測損失は次式で与えられる。


\begin{align*}
R_{\text{err}}(h_{0}) &= \mathbb{E}_{X\sim D_{\mathcal{X}}}   [ \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h_0(X) \neq Y ]  ] ]
\end{align*}

さらに期待値  \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h_{0}(X) \neq Y ]  ]


\begin{align*}
 \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h_{0}(X) \neq Y ]  ] &= \mathrm{Pr}(Y =1 \mid X) \cdot  \mathbf{1} [\beta_X \lt 1/2 ] + \mathrm{Pr}(Y =0\mid X) \cdot \mathbf{1} [\beta_X \geq 1/2 ] \\
&=  \beta_{X} \cdot \mathbf{1} [\beta_X \lt 1/2 ] + (1-\beta_{X})\cdot \mathbf{1} [\beta_X \geq 1/2 ]  \\
&=  \mathrm{min} \{\beta_X, 1-\beta_X \}
\end{align*}

であって、最右辺は  \alpha_{X} の定義そのものである。最終的に


\begin{align*}
R_{\text{err}}(h) \geq R_{\text{err}}(h_{0})
\end{align*}

という不等式を得るが、これは  h_{0} が下限を達成する仮説すなわちベイズ規則であることを示している。

記号をテキストに合わせれば、ベイズ規則は次式で書くこともできる。


\begin{align*}
h_{0} (x) = \mathop{\mathrm{argmax}} \limits_{y \in \mathcal{Y}} \mathrm{Pr}(Y = y \mid X = x)
\end{align*}

さらにまた、


\begin{align*}
\alpha_{X} &= \mathrm{min} \{ \mathrm{Pr}(Y =1 \mid X), \mathrm{Pr}(Y =0 \mid X)\}\\
&=  \mathop{\mathrm{min}} \limits_{y \in \mathcal{Y}} \mathrm{Pr}(Y =y \mid X)\\
&=  1- \mathop{\mathrm{max}} \limits_{y \in \mathcal{Y}} \mathrm{Pr}(Y =y \mid X)
\end{align*}

であるから、テキストの式 (1.2) の0-1損失関数  l_{\text{err}} のもとでのベイズ誤差が次式で与えられることも分かる。


\begin{align*}
R_{\text{err}}^{*} = R_{\text{err}}(h_{0}) &= \mathbb{E}_{X\sim D_{\mathcal{X}}}   [ \mathbb{E}_{Y\sim D_{Y \mid X}}   [ \mathbf{1}  [ h_0(X) \neq Y ]  ] ]\\
&=  \mathbb{E}_{X\sim D_{\mathcal{X}}}   [   \alpha_{X} ] \\
&= 1 -  \mathbb{E}_{X\sim D_{\mathcal{X}}}   \left[  \mathop{\mathrm{max}} \limits_{y \in \mathcal{Y}} \mathrm{Pr}(Y =y \mid X)  \right] 
\end{align*}

参考

統計解析特論-2016(講義資料 4) URL

『統計的学習理論』の第2章のメモ(一様対数の法則, p.32)

一様対数の法則の不等式は以下に示す通り。すなわち、 1-\delta 以上の確率で次式が成り立つ。


\begin{align*}
\mathop{\rm sup} \limits_{g\in \mathcal{G}} \left\{ \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right\}
\leq 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}}
\end{align*}

不等式の左辺における  \text{sup} の中身の符号を反転したものについても同様の不等式が成り立つので、結局  1-\delta 以上の確率で


\begin{align*}
\mathop{\rm sup} \limits_{g\in \mathcal{G}} \left| \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right|
\leq 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (2/\delta)} {2n}}
\end{align*}

が成り立つ。

右辺の  2/\delta は誤植ではない。実際、


\begin{align*}
 A(Z_1, \ldots, Z_n) &= \mathop{\rm sup} \limits_{g\in \mathcal{G}} \left\{ \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right\}\\
 B(Z_1, \ldots, Z_n) &= \mathop{\rm sup} \limits_{g\in \mathcal{G}} \left\{ - \mathbb{E}[g(Z)] + \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right\}\\
 C(Z_1, \ldots, Z_n) &= \mathop{\rm sup} \limits_{g\in \mathcal{G}} \left| \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right|
\end{align*}

とおくと、


\begin{align*}
\text{Pr} \left( A(Z_1, \ldots, Z_n) \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}}  \right) \lt \delta 
\end{align*}

であり、また


\begin{align*}
\text{Pr} \left( B(Z_1, \ldots, Z_n) \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}}  \right) \lt \delta 
\end{align*}

である。和事象の確率に関して成り立つ不等式(劣加法性)を使えば


\begin{align*}
& \text{Pr} \left( C(Z_1, \ldots, Z_n)  \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}} \right) \\
&= \text{Pr} \left( A(Z_1, \ldots, Z_n) \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}} 
\;\lor \right. \\
& \hspace{2cm} \left. B(Z_1, \ldots, Z_n)  \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}}  \right) \\
&\leq  \text{Pr} \left( A(Z_1, \ldots, Z_n)  \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}} \right) \\
& \hspace{5mm} + \text{Pr} \left( B(Z_1, \ldots, Z_n) \gt 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}} \right) \\
& \lt \delta + \delta = 2\delta
\end{align*}

となる。つまり  1-2\delta 以上の確率で


\begin{align*}
\mathop{\rm sup} \limits_{g\in \mathcal{G}} \left| \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right|
\leq 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (1/\delta)} {2n}}
\end{align*}

である。  2\delta \delta に置き換えれば(つまり式の上では  \delta を2で割る)、 1-\delta 以上の確率で


\begin{align*}
\mathop{\rm sup} \limits_{g\in \mathcal{G}} \left| \mathbb{E}[g(Z)] - \frac{1}{n} \sum_{i=1}^{n} g(Z_i) \right|
\leq 2 \mathfrak{R} (\mathcal{G}) + (b-a) \sqrt{\frac{\log (2/\delta)} {2n}}
\end{align*}

が成り立つ。

『統計的学習理論』の誤植

公式の正誤表以外のものを列挙する。

以下の正誤表は私の手元にある第1刷に基づく。

随時追記予定。

p.24

 S_1 S_2の共通部分 →  \textcolor{red}{\text{conv}(S_1)} \textcolor{red}{\text{conv}S_2} の共通部分

p.28


\begin{align*}
\frac{1}{2n} \sum_{y \in \mathcal{Y}} \mathbb{E}_{\sigma} \left[ \mathop{\rm sup} \limits_{f\in \mathcal{G}} 
\sum_{i}(2 \cdot \mathbf{1}[ y = y_i] - 1) f(x_i, y) \right]
\end{align*}


\begin{align*}
\frac{1}{2n} \sum_{y \in \mathcal{Y}} \mathbb{E}_{\sigma} \left[ \mathop{\rm sup} \limits_{f\in \mathcal{G}} 
\sum_{i} \textcolor{red}{\sigma_i} (2 \cdot \mathbf{1}[ y = y_i] - 1) f(x_i, y) \right]
\end{align*}

p.28

最後の不等式 → 最後の等式

p.68 上から9行目

 U: \mathcal{H} \mapsto l(I) →  U: \mathcal{H} \mapsto l^{2}(I)

p.68 上から13行目

 f \in l(I) →  f \in l^{2}(I)