日々のつれづれ

不惑をむかえ戸惑いを隠せない男性の独り言

LARSの整理・その1

LASSOをまとめていると、LARS(最小角回帰法:Least Angle Regression)が出てきたので、これを機会に手順をちゃんと理解しようと思う。

行列と角度、高校や大学で習ったものの、いざ実データで、となるとどうしてもピンと来ない。実装してようやく理解できたので、その記録。

LARSとは

LARSは、LASSOを効率的に計算できるスパースモデリングの一種です。

回帰分析において、変数選択と回帰係数の推定を同時に行うアルゴリズムであり、特徴量が非常に多い高次元データセットに対して効率的に特徴選択できる強みがある。

LARSが生まれた背景

  • LARSはLASSOが持つ課題を開設するために考案されたスパースモデリングの一種。
  • 数値的な最適化手法でパラメータを推定するLASSOがもつ計算コストや収束性の課題に対して、LARSは解析的なパラメータ決定を提供した。
  • 説明変数と残差が成す角度が最も小さい変数を順に選択して、説明変数と残差の角度が等しくなるまで回帰係数の更新を繰り返す。
  • この変数選択法は、LASSOよりも計算コストや収束性などが改善していて、効率的とされている。
  • また、このアプローチは、説明変数間に相関がある場合でも、有効に変数選択や次元削減できる。

LARSの流れ

LARSは最小二乗法をベースにしていて、目的変数との係数がゼロにならない説明変数のセット(アクティブセット)、それを使って求める予測値、この予測値と目的変数の残差を考えます。

次のような流れで計算が進みます。

  1. 初期値として全ての回帰係数を0にする。
  2. 目的変数となす角度が最も小さい(目的変数との相関係数の絶対値が最も大きい)説明変数を求め、この選択した説明変数をアクティブセットにする。
  3. アクティブセットの線形結合で表される方向に沿って、目的変数との残差が小さくなるように回帰係数を更新する。
  4. 更新した後の残差となす角度が最も小さい説明変数を求め、この選択した説明変数をアクティブセットに加える。
  5. アクティブセットの線形結合で表される方向(プロジェクションの方向)に沿って、残差を減らすように回帰係数を更新する。
  6. この残差とアクティブセット以外の列ベクトルとの角度が同じになるまで、4,5を繰り返す。
  7. 残差とアクティブセット以外の列ベクトルとの角度が同じになったら、その列ベクトルをアクティブセットに加える。また、残差ベクトルとアクティブセット内の列ベクトルとの角度が180度になったら、その列ベクトルをアクティブセットから除く。
  8. 全ての説明変数がアクティブセットになるか、残差が0になったら終了する。
\[
\begin{align*}
&\text{Step 1: 初期化} \\
&\quad \beta = \mathbf{0} \quad \text{(全ての回帰係数をゼロベクトルで初期化)} \\
&\quad r = y \quad \text{(残差を目的変数で初期化)} \\
&\quad X_{\text{active}} = \emptyset \quad \text{(アクティブセットを空に初期化)} \\
\\
&\text{Step 2: メインループ} \\
&\quad \text{while } \|\beta\|_0 < p \quad \text{(非ゼロ要素の数が説明変数の総数より小さい間ループ)} \\
\\
&\quad \quad \text{Step 2.1: 最大相関の説明変数を選択} \\
&\quad \quad \text{if } X_{\text{active}} = \emptyset \quad \text{(最初のステップの場合)} \\
&\quad \quad \quad j_{\text{max}} = \text{argmax}_j |\langle X_j, r \rangle| \quad \text{(説明変数\(X_j\)と残差\(r\)の内積の絶対値が最大のインデックス)} \\
&\quad \quad \text{else} \quad \text{(2回目以降のステップの場合)} \\
&\quad \quad \quad X_{\text{active}}^c = X_{\text{active}} \quad \text{(アクティブセットから選択されていない変数のインデックス)} \\
&\quad \quad \quad Z = X^T_{\text{active}} X_{\text{active}}^c \quad \text{(アクティブセットの変数同士の内積行列)} \\
&\quad \quad \quad a = X^T r \quad \text{(各変数と残差の内積)} \\
&\quad \quad \quad j_{\text{max}} = \text{argmax}_j |\langle Z_j, a \rangle| \quad \text{(\(Z_j\)は行列\(Z\)の各列ベクトル)} \\
\\
&\quad \quad \text{Step 2.2: 選択された変数をアクティブセットに追加} \\
&\quad \quad X_{\text{active}} = X_{\text{active}} \cup \{j_{\text{max}}\} \\
\\
&\quad \quad \text{Step 2.3: アクティブセットに対する回帰係数を計算} \\
&\quad \quad \beta_{\text{active}} = (X_{\text{active}}^T X_{\text{active}})^{-1} X_{\text{active}}^T y \\
\\
&\quad \quad \text{Step 2.4: アクティブセットに対する回帰係数を全変数に拡張} \\
&\quad \quad \beta = \mathbf{0} \quad \text{(全ての回帰係数をゼロベクトルで初期化)} \\
&\quad \quad \beta_{\text{active}} \text{の対応する位置に}\beta_{\text{active}} \text{の値をセット} \\
\\
&\quad \quad \text{Step 2.5: 残差の更新} \\
&\quad \quad r = y - X \beta \\
\\
&\quad \quad \text{Step 2.6: 新しい相関を計算} \\
&\quad \quad c = X^T r \\
\\
&\text{Step 3: モデルの終了条件の判定} \\
&\quad \text{if } \|\beta\|_0 = p \text{ または } \|r\|_2 = 0 \quad \text{(全ての説明変数がアクティブセットになるか、残差が0になったら終了)}
\end{align*}
\]

LARSの比較表

LARSはLASSOの解を求めるアルゴリズムなので、RidgeやElastic Netと比較した表はこうなる。

解析法 正則化 変数選択 説明変数間の相関 計算コスト
LARS L1ノルム(LASSO)など 可能(LASSO)など 有利 低い
LASSO L1ノルム 可能 不利 高い
Ridge L2ノルム 不可能 有利 低い
Elastic Net L1ノルム+L2ノルム 可能 有利 高い

LARSは計算コストを抑えてLASSOと同等の性能を出せるらしく、他と比べると使いどころがありそうです。

もちろん、変数選択や次元削減は問題、データの質に強く依存します。 どの方法も万能ではなく、目的に合わせて試行錯誤が必須なことは、言うまでもないですね。

最尤推定の復習・その1

最尤推定法の復習

最尤推定法は、データに確率モデルを当てはめるときに必要になるパラメータを、データから推定する方法です。

私の場合、大学の教養課程で習ったものの在学中に使う機会はなく、就職してからも使う機会なく時間が流れ、40歳を目前としたころ急に必要になりました。習っていても使わぬまま20年過ぎてしまうと、覚えているはずもなく、Rではライブラリで処理できたことから、あいまいな知識のまま、今に至っている方法です。

私にとって最尤推定法を理解できない原因ははっきりしていて、常套句のように出てくる「観測されたデータをもっとも起こりやすい」という表現です。これがしっくりこない。白黒つかない印象がダメです。

とはいえ、これじゃあ、どうしようもないので、今の理解を記録しておこうと思います。

最尤推定法の手順

データ(観測値)を $X_1, X_2, \ldots, X_n$ とすると、尤度関数は $L(\theta) = f(X_1, X_2, \ldots, X_n | \theta)$ と書ける。この式が確率密度関数または確率関数になる。そして、この式の $\theta$ がパラメータで、これを推定するのが最尤推定法です。

また、最尤推定は $L(\theta) = f(X_1, X_2, \ldots, X_n | \theta)$ が成立するときの $\theta$ の最大値を求めるので、 $\hat{\theta} = \arg\max_{\theta} L(\theta)$ となり、この $\hat{\theta}$ を最尤推定値と呼ぶ。

そして、手順おおきく2つに分かれる。

  1. 尤度関数を定義する データに当てはめたい数式が尤度関数になる。尤度関数にはパラメータが必要で、データが連続値なら確率密度関数、離散値なら確率関数のパラメータになる。
  2. 尤度関数を最大化するパラメータの値を決める。 一般的に、最大化には最適化法を使うことが多く、最適化手法には様々な方法がある。

最尤推定の骨格

最尤推定法をRで実装するとき、最適化関数optimがある。最適化方法をオプションで指定できるので、コーディングするときに便利。

で、最尤推定法の手順で考えると、

尤度関数を定義する

パラメータとデータを引数にとるように尤度関数を定義する。

# 尤度関数
likelihood <- function(theta, data) {
  # 使いたい確率密度関数を書き込む
  # パラメータthetaとデータdataを引数にして、尤度を計算する
}

尤度関数を最大化するパラメータの値を決める

最適化方法にはいろいろあって、状況によって決まる。

パラメータの初期値を設定して、反復処理を繰り返して、最適値を探索する。

# 最尤推定
mle <- function(data) {
  # 尤度関数を最大化するパラメータの初期値を決める。
  par0 <- ...

  # 最適なパラメータを探索するための最適化手法を書き込む。
  # パラメータを最適化するためのパラメータ(ハイパーパラメータ)を
  # 定義することもある

  # 最大化した尤度と推定されたパラメータを返す。
  return(list(loglikelihood = ..., estimates = ...))
}

最適化方法はいくつかあるものの、optim関数を使えば、関数の定義が不要になるので、少し楽できる。

注意点としては、optim関数は尤度関数を最小化する関数なので、optim関数に渡す引数(fn)は最小化関数にする必要がある(対数尤度関数)。また、最適化手法やパラメータの初期値などをoptim関数の引数に指定するが、これは試行錯誤で決定する部分が多大にある。

# 最尤推定
mle <- function(data) {
  # optim関数は、デフォルトで尤度関数を最小化するので、
  # 対数尤度関数といった、尤度関数を最小化を考える。
  likelihood <- function(theta) {
    return(-likelihood(theta, data))
  }

  # 最適化関数optimを使ってパラメータを最適化する。
  # 尤度関数(fn)の他、パラメータの初期値(par)や最適化手法(method)も必要。
  result <- optim(par = ..., fn = likelihood, ...)

  # 最小化した負の対数尤度と推定されたパラメータを返す。
  return(list(loglikelihood = -result$value, estimates = result$par))
}

LASSO回帰の復習・その2

Lasso回帰と他の回帰の比較表

Lasso回帰は比較的新しい方法で、L1正則化を行うRidge回帰は少し古く、L1正則化とL2正則化を組み合わせるElastic NetはLasso回帰より後の年代です。 この2つに最小二乗法を加えた2種とLasso回帰を比較します。

  • 最小二乗法
    • 利点: 最も一般的な手法で、解釈が容易で計算が高速。
    • 欠点: 外れ値、多重共線性に弱い。
  • Ridge回帰
    • 誤差の二乗和に重みの二乗和を加えて最小化するL2正則化をで係数を推定する。
    • 利点: 多重共線性に強く、汎化性能が高い。
    • 欠点: 変数選択しないので、解釈性が低い。
  • Lasso回帰
    • 誤差の二乗和に重みの絶対値の和を加えて最小化するL1正則化で係数を推定する。
    • 利点: 疎な解になり、解釈性が高い。
    • 欠点: 相関が高い同士の変数選択が任意になる。
  • Elastic Net回帰
    • 誤差の二乗和に重みの絶対値の和と重みの二乗和の和を加えて最小化するように係数を推定する(L1正則化とL2正則化を組み合わせる)
    • 利点: 疎な解になり、多重共線性に強い。
    • 欠点: ハイパーパラメータのチューニングが必要、L1正則化とL2正則化のバランスが難しい。
手法 発見年代 正則化 特徴選択 パラメータ推定 多重共線性 ノイズ耐性 過学習防止 得意な分類問題
最小二乗法 1800年 なし 非選択 解析的 弱い 弱い なし 線形回帰、多重回帰
Ridge回帰 1970年代 L2正則化 非選択 最小化 強い 強い あり 特徴量の縮小、回帰モデルの安定化
Lasso回帰 1996年 L1正則化 選択 最小化 弱い 強い あり 特徴選択、遺伝子発現データ解析
Elastic Net回帰 2005年 L1正則化 + L2正則化 選択 最小化 強い 強い あり 特徴選択、高次元データの回帰モデリング

LASSO回帰の復習・その1

Lasso(Least Absolute Shrinkage and Selection Operator)回帰とは

Lasso回帰は、統計モデルの予測精度を向上させるための回帰分析手法の一種で、機械学習統計学の分野で広く利用されています。最初の報告は地球物理学で行われましたが、後にRobert Tibshiraniが再発見し、一般化しました。

Lasso回帰の特徴

当時、私は統計モデルの性能が思うように上がらずに悩んでいました。そんなとき、統計学会でLasso回帰の発表を聞き、驚いたのを覚えています。

Lasso回帰はL1正則化という機能を持ち、重要度が低い説明変数の係数(重み)を0にするスパース性を示します。 この重みがパラメータで、重みが0になった説明変数は統計モデルに使われなくなります。このようにパラメータ推定と変数選択が同時に進みます。

しかも、スパース性は変数を少なく抑える効果を持つので、予測モデルの過剰適合(オーバーフィッティング)を抑えることになります。そのため、多重共線性を示す説明変数間の相関が高いケースや、説明変数がサンプル数よりも多いケース、欠損値があるケースなど、今まで使いにくかったデータセットにも適用できる可能性が生まれました。

一方で、Lasso回帰も万能ではなく、パラメータのチューニングが重要、信頼区間の計算が困難という特徴もあります。

とはいえ、L1正則化は強力なので、もとは最小二乗法で定義されていたLasso回帰も、一般化線形モデル、一般化推定方程式、比例ハザードモデル、M推定器などへ拡張され、様々な分野での応用が期待されています。

  1. L1正則化: 重要度が低い説明変数の係数(重み)を0にするスパース性を示します。つまり、Lasso回帰によって推定された重みが0になる変数は、統計モデルにおいて無視されることになります。この性質により、変数の選択とパラメータ推定が同時に行われることが特徴です。

  2. スパース性と変数選択: Lasso回帰はスパースな解を示すため、重要な変数だけが残ります。したがって、Lasso回帰は変数選択の手法としても機能します。これにより、モデルの解釈性が向上し、統計モデルの性能が改善されます。

  3. 過剰適合の抑制: Lasso回帰のスパース性は、変数の数を制限する効果を持ちます。そのため、Lasso回帰は過剰適合(オーバーフィッティング)を抑制する効果があります。特に、多重共線性の問題やサンプル数よりも多い変数のあるデータセット、欠損値のあるデータなど、従来の手法では扱いにくかったデータにも適用できます。

また、Lasso回帰はさまざまな分野で応用されています。例えば、生物学では遺伝子発現データの解析、経済学では株価予測や経済指標の関連性解析、医療では疾患予測やバイオマーカーの選択などです。

分野 応用事例
生物学 遺伝子発現データの解析
経済学 株価予測、経済指標の関連性解析
医療 疾患予測、バイオマーカーの選択
画像処理 特徴選択、物体検出、顔認識
自然言語処理 単語の重要度推定、テキスト分類
プロセス制御 化学反応のモデリング、プロセス最適化
SNS分析 コミュニティ検出、影響力推定
製造業 品質予測、不良要因の特定
環境・エネルギー 気候モデリング再生可能エネルギー予測
保険業 ポートフォリオ最適化、リスク予測
輸送・物流 需要予測、在庫最適化

Lasso回帰の実行手順

Lasso回帰の実行手順は以下の通りです。

  1. データの準備: 解析に適した形式のデータを収集し、特徴量(説明変数)と目的変数に分ける。
  2. データの前処理: 説明変数を標準化や欠損値の処理を行う。
  3. データの分割: データを訓練データとテストデータに分ける。
  4. モデルの構築: Lasso回帰モデルを定義し、正則化パラメーター$\lambda$(ラムダ)を設定する。
  5. モデルの学習: 最適化手法を用いて訓練データからL1正則化項に最適なパラメータを求め、学習済みモデルを得る。
  6. モデルの評価: 平均二乗誤差(MSE)や決定係数($R2$))、AICBICなどを指標に、学習済みモデルの性能を評価する。
  7. 特徴選択と結果解釈: モデルの係数を分析し、重要な特徴の特定や解釈を行う。

回帰や分類に使えるデータセットの一覧

毎回、使うときに探すのも大変だし、使い慣れたものばかりだと正しくコードを書いたのか分かりにくい。

2023年時点の記録を残そうとおもった。

回帰や分類の問題に使えそうなデータセット

Rにあるもの

種類 名前 説明 ライブラリ
回帰 mtcars 自動車の性能に関する情報(燃費など) 組み込み
回帰 airquality 大気中の酸素、二酸化窒素、風速などの情報 組み込み
回帰 swiss スイスの州ごとの特徴(農業、教育、福祉など) 組み込み
回帰 Boston ボストン市の住宅価格に関する情報 MASS
回帰 diamonds ダイヤモンドの特性(カラット、カットなど) ggplot2
分類 iris アヤメの花の品種を分類するための情報 組み込み
分類 BreastCancer 乳がんの診断情報 mlbench
分類 Sonar 音響信号の反射パターンを使用して鉱石か岩かを分類する情報 mlbench
分類 Spam 電子メールがスパムか非スパムかを分類する情報 kernlab
分類 Titanic タイタニック号の乗客の情報と生存結果 Titanic

Rにないもの

種類 名前 説明
回帰 California Housing カリフォルニア州の住宅価格の特徴量とその値
回帰 Energy Efficiency 建物のエネルギー効率に関連する特徴量と目的変数
回帰 Wine Quality 異なる特性に基づいたワインの品質の評価
回帰 Red Wine Quality 赤ワインの特性に基づいた品質の評価
回帰 White Wine Quality 白ワインの特性に基づいた品質の評価
分類 Bank Marketing 銀行の顧客の属性と最終的な預金申し込みの結果
分類 Default Credit Card クレジットカードのデフォルト支払いの予測
分類 Breast Cancer 乳がんの診断情報
分類 Pima Indians Diabetes ピマ・インディアン部族の糖尿病発症の予測
分類 German Credit ドイツの銀行の信用情報と顧客の信用リスク評価

全部を使いこなせていないので、本当に使えるのかはよく分からない

ロジスティック回帰モデルの復習・その3

次に、例題を作って、手順を理解してみる。

5教科の成績から入学試験の合否を予測する

5教科の成績から入学試験の合否を予測するロジスティック回帰モデルを考えてみる。

ここでは、5教科の成績を入力変数、入学試験の合否を出力変数にして、ロジスティック回帰モデルのパラメータを、重み(係数)$w$とバイアス(切片)$b$からなる$\theta$とする。 そして、入力変数を反復学習して、ロジットを計算し、シグモイド関数を適用して確率を求め、コスト関数からパラメータを推定する。

データセットの準備

20人分のデータセットを作成して、最初の15人分をロジスティック回帰モデルの構築に使う学習セット、最後の5人分を評価用いる新しいデータセットにする。

このとき、5教科の成績を0から100点からなる入力変数、入学試験の合否を合格なら1、不合格なら0にした出力変数にする。なお、実際には、データのスケーリングや正規化などの前処理も必要ですが、ここでは省略する。

データセット 合否
1 学習 75 82 90 68 78 1
2 学習 92 88 77 80 85 1
3 学習 63 70 75 58 65 0
4 学習 81 78 86 92 88 1
5 学習 70 65 72 68 74 0
6 学習 95 92 88 97 90 1
7 学習 72 75 68 64 70 0
8 学習 88 92 84 78 80 1
9 学習 60 58 63 52 55 0
10 学習 78 80 76 82 85 1
11 学習 68 72 64 58 62 0
12 学習 83 87 92 89 94 1
13 学習 76 70 72 78 80 1
14 学習 55 62 58 50 60 0
15 学習 90 85 92 88 94 1
16 評価 72 74 68 76 80 1
17 評価 62 58 65 60 64 0
18 評価 80 82 85 78 84 1
19 評価 67 72 66 58 70 0
20 評価 94 90 88 92 96 1
X = c(
  75, 82, 90, 68, 78, 1,
  92, 88, 77, 80, 85, 1,
  63, 70, 75, 58, 65, 0,
  81, 78, 86, 92, 88, 1,
  70, 65, 72, 68, 74, 0,
  95, 92, 88, 97, 90, 1,
  72, 75, 68, 64, 70, 0,
  88, 92, 84, 78, 80, 1,
  60, 58, 63, 52, 55, 0,
  78, 80, 76, 82, 85, 1,
  68, 72, 64, 58, 62, 0,
  83, 87, 92, 89, 94, 1,
  76, 70, 72, 78, 80, 1,
  55, 62, 58, 50, 60, 0,
  90, 85, 92, 88, 94, 1,
  72, 74, 68, 76, 80, 1,
  62, 58, 65, 60, 64, 0,
  80, 82, 85, 78, 84, 1,
  94, 90, 88, 92, 96, 1
)
X <- data.frame(matrix(X, ncol = 6, byrow = TRUE))
colnames(X) = c("英", "数", "国", "理", "社", "合否")

y <- factor(X$"合否")
X <- X[, 1:5]
new_X <- X[-c(1:15), ]
new_y <- y[-c(1:15)]
X <- X[1:15, ]
y <- y[1:15]

ロジスティック回帰モデルの考え方

回帰を考えるとき、予測対象が合格(1)、不合格(0)1と0の二値だと扱いにくい。予測対象を合格確率と考えて、テストの点数が高いほど1に近づき、低いほど0に近づくと考えます。

線形予測子を作る

まず、入力変数$x$を5教科の点数、$P(x)$を合格確率と考えます。

このとき、5教科の点数を個別に見るよりも、一つにまとめて、その大小で合格確率を考える方が楽です。そこで、5教科の点数をまとめる式$z$を考えます。

式$z$は係数が$w$、切片が$b$になる多項式で、出力変数$y$の予測に用いる線形予測子になります。

$$ z = w^{T}x + b $$

ここで、$w$と$b$は、ロジスティック回帰モデルのパラメータ$\theta$の重みとバイアスになります。

ロジット関数を作る

次に、出力変数$y$は合格確率なので範囲に上下限があります。これは使いにくいので、オッズ比の対数をとるロジット関数を考えます。

$$ logit(y)=log \frac{y}{1-y} $$

$logit(y)$の予測にすることで、予測範囲は$- \infty \sim \infty$になります。

リンク関数を設定する

次に、多項式$z$が出力変数$y$のロジット関数$logit(y)$になると考えます。このとき、ロジット関数は$x$と$y$の関係を示すリンク関数になります。

$$ logit(y) = z $$

これを変換すると、ロジスティック関数(シグモイド関数)になります。

$$ y = \frac{1}{1+e^{(w^{T}x + b)}} $$

パラメータ推定、学習、誤差計算

実計算では、パラメータ$\theta$を繰り返し修正しながら、作成したロジスティック回帰モデルの予測値と実測値との乖離を最小化して行きます。

このパラメータ$\theta$の修正は、最尤推定法でパラメータ$\theta$を推定する、予測値と実測値のクロスエントロピー誤差を求める、勾配降下法でこの誤差を最小化する流れで進みます。

# ロジット関数の定義
logits <- function(X, w, b) X %*% w + b

# シグモイド関数の定義
sigmoid <- function(z) 1 / (1 + exp(-z))

# クロスエントロピー誤差(コスト関数)の定義
cost <- function(y, y_pred) {
  sum(y * log(y_pred) + (1 - y) * log(1 - y_pred))
}

この定義では、sigmoid(logit(X, w, b))が予測値の計算になります。

また、クロスエントロピー誤差のコスト関数はyが実際の合否、y_predが合否の予測値で、y_predyに近いほど小さくなります。まず、y * log(y_pred)で合格(1)の予測確率の対数を求め、次に、(1 - y) * log(1 - y_pred)で不合格(0)の予測確率の対数を求めます。次に、この2つを足し合わせ、最後に合計して全体誤差を求めています。

# パラメータの初期化
w <- matrix(0, nrow = ncol(X), ncol = 1)
b = matrix(0, nrow = 1, ncol = 1)

# ハイパーパラメータの設定
learning_rate <- 0.01
num_iterations <- 100

# 最尤推定法によるパラメータの最適化
cost_iter = c()
for (i in 1:num_iterations) {
  # 予測値の計算
  y_pred <- sigmoid(logit(X, w, b))

  # 勾配の計算
  dw <- (1 / length(x)) * sum((y_pred - y) * x)
  db <- (1 / length(x)) * sum(y_pred - y)

  # パラメータの更新
  w <- w + learning_rate * dw
  b <- b + learning_rate * db

  cost_iter = c(cost_iter, cost(y, y_pred))
}

ここでは、勾配降下法の反復計算でパラメータ$\theta$を最適化しています。

まず、予測値y_predを求めて、実際の合否yとの差分から勾配を計算します。次に、この勾配に学習率learning_rateを掛け合わせて、パラメータを更新する程度を決めます。最後に、これを重み$w$とバイアス$b$に加えて、パラメータを更新します。指定の回数num_iterationsの反復を行って、最後に得た重み$w$とバイアス$b$が最適化したパラメータ$\theta$と考えます。

このとき、反復計算の各回のクロスエントロピー誤差をcost_iterに貯めています。これは、計算後に誤差の値の減少を確認するためです。この値が減少していなければ、学習率や反復回数の見直しが必要です。

また、誤差の減少が重要であれば、固定値でなく、誤差が指定の値以下になるまで反復計算する方法もあります。ですが、過度な反復計算は過学習(overfitting)のリスクが上げ、未知のデータに対する汎化性能を下げてしまいます。ここには注意が必要です。

一般的には、反復計算の回数を適切な設定の他、正則化(regularization)や交差検証(cross-validation)などの手法を使うことになるでしょう。

予測

学習が終了した$\theta$を使用して新しいデータを予測します。ロジスティック回帰モデルは確率を出力するため、0.5以上なら1と判定し、0.5未満なら0と判定します。

# 正解率、適合率、再現率、F1スコアを計算
pred_score = function(predictions, true_labels){
  # 正解率(Accuracy)の計算
  accuracy <- sum(predictions == true_labels) / length(predictions)

  # 適合率(Precision)と再現率(Recall)の計算
  TP <- sum(predictions == 1 & true_labels == 1)  # True Positive
  FP <- sum(predictions == 1 & true_labels == 0)  # False Positive
  FN <- sum(predictions == 0 & true_labels == 1)  # False Negative
  precision <- TP / (TP + FP)
  recall <- TP / (TP + FN)

  # F1スコア(F1 Score)の計算
  f1_score <- 2 * (precision * recall) / (precision + recall)

  score = list(accuracy, precision, recall, f1_score)
  return(score)
}

# 予測結果と真の値のサンプルデータ(0: Negative, 1: Positive)
pred_y = sigmoid(logit(new_X, w, b))
pred_y <- ifelse(pred_y >= .5, 1, 0)

pred_score(prediction=pred_y, true_labels=new_y)

データを理解する・その3「特異値分解、固有値」

これも、2020年にまとめていた記事。これが最後 放置していたままも嫌なので、公開します。

ここから

長い年末休暇を利用して、理解した記録です。

特異値分解

特異値分解は、行列を 3 つの行列の積に分解する。

$$ X = USV^{T} $$

このとき、$U$を右得意行列、$V$を左得意行列、$U$を特異値行列とよぶ。

$U$と$V$は正方行列で、$U$が$[n \times n]$、$V$が$[m \times m]$のとき、$V$は$[n \times m]$の行列になる。

$S$の要素は次で$s_{1 \dots k}$が特異値になる。また、$SS^{T} = S^{T}S = S^{2}$の性質を持つ。


S = \left(
  \begin{array}{cccc}
  s_{1}  & 0      & \dots  & 0     \\
  0      & s_{2}  &        & 0     \\
  \vdots &        & \ddots & \vdots\\
  0      & 0      &        & s_{k} \\
  0      & 0      & \dots  & 0     \\
  \vdots & \vdots & \vdots & \vdots\\
  \end{array}
\right)

固有値

固有値は、特異値分解の$U$と$V$で$V^{T} = U$となるときに成立する。

$$ A = P \Lambda P^{T} $$

このときの$\Lambda$は次で$\lambda_{1 \dots k}$が固有値


\Lambda = \left(
  \begin{array}{cccc}
  \lambda_1 & 0         & \dots  & 0         \\
  0         & \lambda_2 &        & 0         \\
  \vdots    &           & \ddots & \vdots    \\
  0         & 0         &        & \lambda_k \\
  \end{array}
\right)