「川越達也の抜き打ち発掘レストラン」採点結果の平均点を計算
テレビ番組お願い!ランキングで不定期に川越達也の抜き打ち発掘レストラン！という企画が放送されています。このページの内容紹介によると
美食アカデミーでもお馴染みの川越シェフと
進行のハライチが、
街の隠れた名店を発掘するべく突撃取材！
その場で直接取材交渉し、撮影OKの場合、
その店に自慢の料理を出してもらい川越シェフが試食！
100点満点でガチ採点しランキングを作成する企画！
とのことで、レストランで出された料理を採点していくという番組です。発掘場所を変えながらこれまで何度か放送されており、１放送あたりおよそ４店のレストランの料理を採点しています。
ある程度の採点データがあるので平均値を算出してみました。
採点結果の平均値
まず採点結果の分布を見ます。
９０点、９５点、１００点のキリが良い数字が多いことから、厳密な採点基準があるのではなく川越シェフによるフィーリングで採点されているようです。
単純に計算すれば平均点は９３．８点となります。
しかし、発掘場所や取材対象レストラン、川越シェフの気分などの不確定で偶然的な要素も考慮して平均を求めた方が良いので、ここではブートストラップ法で平均値を求めてみます。
実際に計算すると平均値の９５％信頼区間（BCa 法 (bias-corrected and accelerated percentile method)による）は、92.79～94.87点となりました。
したがって、川越シェフの採点結果から、統計学に基づいて次のように料理の味を判断できます（あくまでも川越シェフの主観ですが）。
- ９２点以下の料理は平均より美味しくないといえます
- ９３点から９４点の料理は平均的な味といえます
- ９５点以上の料理は平均より美味しいといえます
ところで、上の図をみると８０点のデータだけ他の得点と比べて低すぎるので、外れ値（異常値）の可能性があります。
そこで８０点のデータだけを無視して上と同じ計算をしました。
単純な平均点は９４．０点でした。
ブートストラップ法による平均点の信頼区間は、93.03～95.04点となりました。
したがって
- ９３点以下の料理は平均より美味しくないといえます
- ９４点から９５点の料理は平均的な味といえます
- ９６点以上の料理は平均より美味しいといえます
といえます。
８０点を外れ値と見るかどうかで結果が若干異なるのですが、９６点以上を獲得した料理であれば平均よりは美味しいと自信をもって言えそうです。
参考サイト
得点データは次のサイトから取得しました。
サポートベクターマシンとは［カーネル法による非線形サポートベクターマシン］
ここからはこれまで述べてきたサポートベクターマシンにカーネル法を適用することにより非線形サポートベクターマシンへ拡張することを考えます。
カーネル法の導入
これまで述べてきたサポートベクターマシーン分離面が超平面であることを前提としていました。しかし、実際の問題では正例データと負例データの境界が超平面であるよりは、複雑に入り組んだ超曲面である可能性が高いことが想定されます。
このようなデータに、これまで述べてきたようなクラス境界を超平面とするサポートベクターマシーンを適用しても、高い分類性能を期待することは難しそうです。
たとえば下図のような単純なケースでさえ正例データ（○）と負例データ（×）を分ける直線は存在しないため、１００％の分類性能は達成できません。
しかし、このようなデータでも線形分離が可能になるような別の空間へ変換できれば、変換先の空間ではクラス境界が超平面になるのでサポートベクターマシーンを用いることが可能になります（下図参照）。元のデータが存在する空間のことを入力空間、変換先の空間のことを特徴空間といいます。一般的に入力空間の次元$n$より特徴空間の次元$M$の方が大きくなります（同じ場合もあります）。実際にこのような写像の例をあとで示します。
入力空間のデータ$x$から特徴空間のデータ$z$へ写像する関数$\phi(x)$とすると$z$は
$$
\begin{align}
z&=(z_1,\ldots,z_M)^T\\
&=(\phi_1(x),\ldots,\phi_M(x))^T
\end{align}
$$
と表現できますので、変換後のデータ$z$に対してソフトマージンサポートベクターマシンを適用します。
特徴空間の例
ここで一旦話を中断し、入力空間では線形分離できないデータでも高次元の特徴空間に写像すると線形分離可能になる例を示しておきます。
下図のような２次元訓練データがあるとします。赤データを正例、緑データは負例としておきます。このデータを正例と負例に分割する直線が存在しないことは明らかです。しかし曲線であれば両クラスを分離する事ができます。図の黒い円（$x_1^2+x_2^2=1$）がその一例です。
ここで、２次元入力空間の各点$x=(x_1, x_2)$を３次元特徴空間の点$z=(z_1, z_2, z_3)$に写像する変換
$$
z=\phi(x_1, x_2)=(x_1^2, x_2^2, \sqrt{2}x_1x_2)
$$
を考えます。この写像$\phi$で図の黒い円上のデータを変換すると、３次元特徴空間の平面$w^Tz+b=0$の上に乗ることがわかります。実際$w=(1, 1, 0)^T, b=-1$とおいた場合の平面の式は$z_1+z_2+0z_3-1=x_1^2+x_2^2-1=0$となるからです。
実際にこの図のデータを関数$\phi$で３次元特徴空間に写像すると、下図のように赤データと緑データを平面で分離できるようになりました。この特徴空間でサポートベクターマシンを適用すれば１００％の分類精度が得られます。
YouTubeに入力空間から特徴空間への写像を説明した動画がありましたので紹介します。$\phi(x)$の関数形は上の例と異なるようですが、イメージ的には同じです。
特徴空間での定式化
話を戻して特徴空間でのソフトマージンサポートベクターマシンを構成しましょう（$x$を$z$に置き換えるだけです）。
特徴空間でのソフトマージンサポートベクターマシン
訓練データ
$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}$
と、特徴空間への写像$z=\phi(x)=(\phi_1(x),\ldots,\phi_M(x))$が与えられたとき、
$$
\begin{align}
\mathcal{F}&= \{(z_1, y_1), \ldots, (z_l, y_l)\}, (z_i, y_i) \in R^M \times \{-1, 1\},\\
z_i&=(\phi_1(x_i),\ldots,\phi_M(x_i))^T, i=1,\ldots,l
\end{align}
$$
は線形分離可能とします。次の最適化問題を解き、$\alpha$を求めす。
- 目的関数
$$
\sum_{i=1}^{k}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jz_i^Tz_j
$$
を制約条件のもとで$\alpha$に関して最小化制約条件：
$$
C \geqq \alpha_i\geqq 0, i=1,\ldots,k\\
\sum_{i=1}^{l}\alpha_iy_i=0
$$決定関数は
$$
\hat{f}(z) = \text{sgn}(\sum_{\mathcal{V}}\alpha_iy_iz_i^Tz + b)
$$
となります。ここで$b$は
$$
\begin{align}
b&=-\frac{1}{2}\sum_{\mathcal{V}}(\alpha_iy_iz_i^Tz_{+}+\alpha_iy_iz_i^Tz_{-})
\end{align}
$$
$\mathcal{V}$はサポートベクターの集合です。
さて、一般に特徴空間の次元$M$は入力空間の次元より高くなります。場合によっては$M=\infty$になることもあります。高次元になるほど内積$z^T_iz_j$の計算量は増加し、関連データを保持するためのメモリも増加していきます。
ためしに多項式写像によって特徴空間に変換されたデータを保持するために必要とされるメモリ量を帰納的に見積もってみます。
$$
\begin{align}
M&={}_dC_2\\
&=\frac{d!}{(d-2)!2!}\\
&=\frac{d(d-1)}{2}
\end{align}
$$
となります。
たとえば入力空間が１００次元（$d=100$）であれば、特徴空間の次元は約５０００次元($100\times(100-1)/2$)になります。入力データがテキストデータのときには入力空間の次元が数千～数万になることもあります。もし１万次元のデータを写像すると特徴空間の次元は５千万次元となります。
サポートベクターマシンを適用するまえに入力空間のデータを特徴空間に変換しておくと、訓練データ１件を保持するのに必要なメモリ量は約４００MB(50,000,000*8=400,000,000)になります。訓練データが１件ということはないので、あっという間にメモリが不足してしまいます。
このように特徴空間に明にデータを写像するような状況を回避するのがこれから述べるカーネル法です。
前述したサポートベクターマシンのアルゴリズムをよく見ると特徴空間でのデータ$z$はすべて内積として、すなわち$z_i^Tz_j=\phi^T(x_i)\phi(x_j)$として現れています。
実はカーネルと呼ばれる関数$K$を使えば、特徴空間での内積計算$\phi(x)^T\phi(y)$を回避できます。すなわち
$$
K(x,y)=\phi(x)^T\phi(y)
$$
とできます。ただ無条件にこのようなカーネルを構成できるのではなく、マーサー(Mercer)の条件を満たす必要があります（詳細は参考文献参照）。
また、カーネルを用いるメリットは特徴空間への写像$\phi(x)$の具体的な形がわからなくても内積が計算できます。さらに、高次元の特徴空間で内積を計算するのではなく、入力空間でカーネルを計算すればよいので計算量とメモリを大幅に削減できます。
特徴空間におけるサポートベクターマシンのアルゴリズムにあらわれる内積をすべてカーネルで置き換えてしまえば内積計算を回避できます。これをカーネルトリックと呼びます。カーネルトリックにより特徴空間での変数を陽に扱わない手法のことをカーネル法と呼びます。
訓練データ
$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}$
が与えられたとき、次の最適化問題を解き$\alpha$を求めます。
- 目的関数
$$
\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jK(x_i,x_j)
$$
を制約条件のもとで$\alpha$に関して最小化
- 制約条件
$$
C \geqq \alpha_i\geqq 0, i=1,\ldots,l\\
\sum_{i=1}^{l}\alpha_iy_i=0
$$決定関数は
$$
\hat{f}(x) = \text{sgn}(\sum_{\mathcal{V}}\alpha_iy_iK(x_i,x) + b)
$$
となります。ここで$b$は
$$
\begin{align}
b&=-\frac{1}{2}\sum_{\mathcal{V}}(\alpha_iy_iK(x_i, x_{+})+\alpha_iy_iK(x_i,x_{-}))
\end{align}
$$
$\mathcal{V}$はサポートベクターの集合です。
これで分類問題に対するサポートベクターマシンの最終形が完成しましました。
カーネルの例
サポートベクターマシンの特徴の一つは、解きたい問題に応じてカーネルを適切に選び汎化能力を向上できることです。実際、問題の種類に応じて様々なカーネルが開発されています。
ここでよく用いられるカーネル$K(x,x')$をいくつか紹介します。
- 線形カーネル
もとの訓練データが線形分離可能であればわざわざ特徴空間に写像する必要はないので、このような場合には次の線形カーネルを用いることにより実質的にはカーネル法を用いない線形のサポートベクターマシンになります。
$$
K(x,x')=x^Tx'
$$
テキストデータのように大規模で疎なデータ（ほとんどの列が０であるようなデータ）でよく用いられます。
- 多項式カーネル
次数を$d$とする多項式カーネルは次式で与えられます。パラメータ$d, scale, offset$はユーザーが指定する必要があります。
$$
K(x,x')=(scale \times x^Tx'+offset)^d
$$
画像を分類するときによく用いられます。なお$offset$は通常１とします。$d=1, scale=1, offset=0$とすれば線形カーネルとなります。
- ラジアル基底関数カーネル（RBFカーネル、ガウシアンカーネル）
ラジアル基底関数(RBF, radial basis function)カーネル（ガウシアンカーネルとも呼ばれます）は次式で与えられます。パラメータ$\sigma(>0)$はユーザーが指定します。
$$
K(x,x')=\exp(-\sigma||x-x'||^2)
$$
データに関する事前知識がない場合に用いられる汎用的なカーネルです。最もよく用いられるカーネルの一つです。
- 双曲線正接カーネル
双曲線正接カーネル(hyperbolic tangent kernel)は次式で与えられます。パラメータ$scale, offset$はユーザーが指定します。
$$
K(x,x')=\tanh(scale \times x^Tx'+offset)
$$
おもにニューラルネットワークの代わりとして用いられます。
この他にも文字列カーネルなど色々ありますが省略します。
最後に$K(x,y)=\phi(x)^T\phi(y)$となるカーネル$K$を実際に構成できることを示す例を挙げておきます。
つぎの写像を考えます。
$$
z=\phi(x_1, x_2)=(x_1^2, x_2^2, \sqrt{2}x_1x_2)
$$
特徴空間での内積$\phi(x)^T\phi(y)$を計算してみると
$$
\begin{align}
\phi(x)^T\phi(y)&=(x_1^2, x_2^2, \sqrt{2}x_1x_2)^T(y_1^2, y_2^2, \sqrt{2}y_1y_2)\\
&=x_1^2y_1^2+x_2^2y_2^2+2x_1x_2y_1y_2\\
&=(x_1y_1+x_2y_2)\times(x_1y_1+x_2y_2)\\
&=(x_1y_1+x_2y_2)^2\\
&=(x^Ty)^2\\
&=(1 \times x^Ty+0)^2
\end{align}
$$
となり、これは多項式カーネルで次数$d=2$、$scale=1$、$offset=0$とおいた式と等しいことがわかります。
この例のように特徴空間で内積$\phi(x)^T\phi(y)$をそのまま計算する代わりに、入力空間で$(x_1y_1+x_2y_2)^2$を計算すれば良いことがわかります。
ほかにも色々と書きたいと思っていたのですが、案外時間がかかりそうなので、このへんで終わります。
参考文献
- Nello Cristianini、John Shawe‐Taylor、大北 剛 (翻訳) 、サポートベクターマシン入門、 共立出版
- 小野田 崇、サポートベクターマシン (知の科学) 、オーム社
サポートベクターマシンとは［ソフトマージンサポートベクターマシン］
スラック変数の導入
スラック変数を導入すると、訓練データの各データが支持超平面から分類超平面のほうにどの程度はみ出したかを測ることができます。別の表現をすれば、はみ出したデータを無視して支持超平面を構成した結果として発生する誤差の程度を測る変数です。
スラック変数を$\xi=(\xi_1, \ldots, \xi_l)^T$と表記すると、スラック変数は誤差の大きさを表しているので
$$
\xi_i\geqq 0, i=1,\ldots, l
$$
です。
下図のような状況を考えます。このデータの配置では正例データと負例データを分ける超平面はもはや存在しないことがわかります。ピンク色の点は誤差データと考え、それ以外のデータから支持超平面を構成しています。誤差データに対してはスラック変数$\xi_i$でその大きさを考慮しています。
支持超平面からはみ出したデータに対応したスラック変数のみ$\xi_i> 0$となり、それ以外は$\xi_i= 0$となることに注意してください。
ハードマージンサポートベクターマシンで出てきた最適化問題にスラック変数を導入すると次のようになります。
スラック変数を導入した最適化問題（主問題）
訓練データ
$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}$
と定数$C(>0)$が与えられたとき、制約条件のもとで次の目的関数を最小化します。
- 目的関数
$$
\frac{||w||^2}{2}+C\sum_{i=1}^{l}\xi_i
$$
- 制約条件
$$
\begin{align}
y_i(w^Tx_i+b) &\geqq 1-\xi_i, i=1,\ldots, l\\
\xi_i &\geqq 0, i=1,\ldots, l
\end{align}
$$
この最適化問題についていくつかコメントしておきます。
- ハードマージンサポートベクターマシンでの最適化問題との違いについて
目的関数には誤差項($\sum_{i=1}^{l}\xi_i$)が新たに加わり、制約条件にはスラック変数の制約が加わっています。
ここで誤差項の意味を考えてみます。
もし目的関数に誤差項を加えない（$C=0$）と、制約条件を満たしつつ、マージンをいくらでも大きくできてしまいます。たとえば、訓練データ全てを誤差データとみなしまうことによりマージンを最大化できます。このような好ましくない状況を避けるためにこの誤差項がペナルティーとして必要であることが分かります。
ちなみに、この最適化問題では誤差項は$\xi$の一次式ですが、$p$次式を用いて次のように目的関数を構成することも可能です。
$$
\frac{||w||^2}{2}+\frac{C}{p}\sum_{i=1}^{l}\xi_i^p
$$
- 定数Cについて
$C(>0)$は誤差項($\sum_{i=1}^{l}\xi_i$)と$||w||^2/2$の大きさのバランスまたはトレードオフを調整するパラメータでユーザーが指定する必要があります*1。
ここで、$C$の効果を考えてみます。
まず、$||w||^2/2$と$\sum_{i=1}^{l}\xi_i$の大小関係を見てみます。
$||w||^2$を小さく（すなわちマージンを大きく）すると支持超平面の反対側にはみ出すデータが多くなります。それに伴い誤差項は大きくなります。逆に$||w||^2$を大きく（すなわちマージンを小さく）すると支持超平面の反対側にはみ出すデータが少なくなるため誤差項は小さくなります。このように目的関数の一方の項を小さくしようとすればもう片方の項が大きくなるので、そのトレードオフをパラメータ$C$で調整します。
つぎに、$C$と$\sum_{i=1}^{l}\xi_i$の関係を見てみます。
最適化問題の目的は目的関数の値を小さくすることなので、１項目だけでなく２項目の$C\sum_{i=1}^{l}\xi_i$も小さくする必要があります。
もし$C$に大きな値を設定すると$C\sum_{i=1}^{l}\xi_i$は当然大きくなりますが、その程度を緩和するためには、$\sum_{i=1}^{l}\xi_i$を小さくする必要があります。$\xi_i \geqq 0$なので、これは$\xi_i$を全体的に小さくするすることになります。すなわち、$C$を大きくすると、誤差データをあまり許容できないためマージンの幅が狭くなります。
逆に$C$に小さい値を指定すると$\sum_{i=1}^{l}\xi_i$をある程度大きくしても$C\sum_{i=1}^{l}\xi_i$はそれほど大きくなりません。すなわち、マージン幅を広く取り誤差データを許容できる程度が大きくなります。
このようにユーザー設定パラメータである$C$の値に応じて分類超平面の位置が変わるため、最も高い分類性能が得られる$C$の値を探すことになります。
ところで、パラメータに文字$C$が用いられることが多いため、ソフトマージンサポートベクターマシンのことをC-SVC(C-support vector classification)と呼ぶこともあります*2。
- 決定関数について
スラック変数を導入しても決定関数はそのままで
$$
\hat{f}(x) = \text{sgn}(w^T x + b)
$$
です。この判別式に、分類超平面より負例（正例）側にある正例（負例）データを代入すると負例（正例）データと判別されるため、このデータは誤分類されます。しかし、もともとこのようなデータを誤差データとみなしていたので、誤分類されたとしても影響は小さいと考えられます。
最適化問題の解法
ハードマージンサポートベクターマシンの場合と同様に主問題の最適化問題を双対問題へ変換しましょう。
ソフトマージンサポートベクターマシンの場合のラグランジュ関数は次で定義されます。
$$
L(w, b, \xi, \alpha, \beta)=\frac{||w||^2}{2}+C\sum_{i=1}^{l}\xi_i-\sum_{i=1}^{l}\alpha_i(y_i(w^Tx_i+b)-1+\xi_i)-\sum_{i=1}^{l}\beta_i\xi_i
$$
ここで、$\alpha=(\alpha_1,\ldots, \alpha_k)^T, \beta=(\beta_1,\ldots, \beta_m)^T$はラグランジュ乗数です。
この場合のクーン・タッカーの定理の関係式はつぎのように表現されます。
$$
\begin{align}
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial w}&=0,\\
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial b}&=0,\\
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial \xi}&=0,\\
\alpha_i(y_i(w^Tx_i+b) -1+\xi_i)&=0, i=1,\ldots,l,\\
\beta_i \xi_i&=0, i=1,\ldots,l,\\
\alpha_i&\geqq 0, i=1,\ldots,l,\\
\beta_i&\geqq 0, i=1,\ldots,l,\\
\xi_i&\geqq 0, i=1,\ldots,l,
\end{align}
$$
４番目と５番目の式がKKT相補性条件となります。
まずラグランジュ関数の偏微分を計算します。
$$
\begin{align}
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial w}&=w-\sum_{i=1}^{l}\alpha_iy_ix_i\\
&=0\\
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial b}&=-\sum_{i=1}^{l}\alpha_iy_i\\
&=0\\
\frac{\partial L(w, b, \xi, \alpha, \beta)}{\partial \xi_i}&=C-\alpha_i-\beta_i \\
&=0
\end{align}
$$
整理すると
$$
\begin{align}
w&=\sum_{i=1}^{l}\alpha_iy_ix_i,\\
\sum_{i=1}^{l}\alpha_iy_i&=0,\\
C&=\alpha_i+\beta_i, i=1,\ldots,l\\
\end{align}
$$
これらの式をラグランジュ関数に代入すれば主問題に対する双対問題の目的関数が得られます。
$$
\begin{align}
L(w, b, \xi, \alpha, \beta)&=\frac{w^Tw}{2}+C\sum_{i=1}^{l}\xi_i-\sum_{i=1}^{l}\alpha_i(y_i(w^Tx_i+b)-1+\xi_i)-\sum_{i=1}^{l}\beta_i\xi_i\\
&=\frac{(\sum_{i=1}^{l}\alpha_iy_ix_i)^T(\sum_{i=1}^{l}\alpha_iy_ix_i)}{2}+\sum_{i=1}^{l}(\alpha_i+\beta_i)\xi_i\\
&-\sum_{i=1}^{l}\alpha_i\{y_i((\sum_{i=1}^{l}\alpha_iy_ix_i)^Tx_i+b)-1+\xi_i\}-\sum_{i=1}^{l}\beta_i\xi_i\\
&=\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j+\sum_{i=1}^{l}\alpha_i\xi_i+\sum_{i=1}^{l}\beta_i\xi_i\\
&-\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j-b(\sum_{i=1}^{l}\alpha_iy_i)+\sum_{i=1}^{l}\alpha_i-\sum_{i=1}^{l}\alpha_i\xi_i-\sum_{i=1}^{l}\beta_i\xi_i\\
&=\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j (\because \sum_{i=1}^{l}\alpha_iy_i=0)
\end{align}
$$
また$C=\alpha_i+\beta_i$を変形すると
$$
\beta_i=C-\alpha_i, i=1,\ldots,l
$$
となり、ラグランジュ乗数$\beta_i \geqq 0$であることから、最終的に$C \geqq \alpha_i \geqq 0, i=1,\ldots,l$が得られます。
以上をまとめると$\alpha$に関する次の最適化問題（双対問題）が得られます。
最適化問題（双対問題）
- 目的関数
$$
\sum_{i=1}^{k}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j
$$
を次の制約条件のもとで$\alpha$に関して最小化
- 制約条件
$$
C \geqq \alpha_i\geqq 0, i=1,\ldots,l\\
\sum_{i=1}^{l}\alpha_iy_i=0
$$
この最適化問題を解いて$\alpha$を求めれば、$w$を計算できます。
最適化問題（主問題）でスラック変数$\xi$を導入したにもかかわらず、この双対問題にはスラック変数が現れないことに注意しましょう。そのため、目的関数の式はハードマージンサポートベクターマシンのときと同じになっています。
ただしハードマージンサポートベクターマシンとの相違点は、制約条件$C \geqq \alpha_i\geqq 0, i=1,\ldots,l$にパラメータ$C$が現れていることです（ソフトマージンでは$\alpha_i\geqq 0, i=1,\ldots,l$でした）。
この制約条件の式$C \geqq \alpha_i\geqq 0$から$\alpha_i$の取りうる値に関して次の３つの場合が考えられます。
- $\alpha_i=0$のとき
$C=\alpha_i+\beta_i$なので$C=\beta_i$となり、$C>0$であることに注意すると$\beta_i>0$となります。このときKKT相補条件の式$\beta_i \xi_i=0$から$\xi_i=0$となります。すなわち双対問題の解$\alpha_i=0$となるデータでは$\xi_i=0$なので、訓練データ$x_i$は支持超平面の正例側または負例側にある（支持超平面からはみ出していない）ことが分かります。したがってこのデータは正しく分類されます。
- $0<\alpha_i < C$のとき
$\beta_i=C-\alpha_i$なので$\beta_i>0$となります。このとき、KKT相補条件の式$\alpha_i(y_i(w^Tx_i+b) -1+\xi_i)=0$、$\beta_i \xi_i=0$から、それぞれ$y_i(w^Tx_i+b) -1+\xi_i=0$、$\xi_i=0$となるので、両式から$y_i(w^Tx_i+b) -1=0$が成り立つことが分かります。これは訓練データ$x_i$が支持超平面にちょうど乗っているので、$x_i$はサポートベクターとなります。当然ながらこのデータは正しく分類されます。
- $\alpha_i = C$のとき
$C=\alpha_i+\beta_i$なので$\beta=0$となります。このとき、KKT相補条件の式$\alpha_i(y_i(w^Tx_i+b) -1+\xi_i)=0$、$\beta_i \xi_i=0$から、それぞれ$y_i(w^Tx_i+b) -1+\xi_i=0$、$\xi_i \geqq 0$が得られます。すなわち、訓練データ$x_i$は支持超平面の反対側にはみ出していることになります。このようにはみ出しているデータもサポートベクターと呼ばれます。
$0 \leqq \xi_i \leqq 1$のときは、$x_i$は支持超平面から反対側にはみ出していますが分類超平面は超えていない（下図の点$x_A$）ため正しく分類されます。しかし$\xi_i \geqq 1$になると、$x_i$は支持超平面だけではく分類超平面も超えてはみ出している（図の点$x_B$、$x_C$）ためこのデータは誤分類されることになります。
$w$の計算量削減
ハードマージンサポートベクターマシンの場合と同様に、$\alpha_i=0$となるデータはサポートベクターではないことを利用すれば、$w$の計算量を削減できます。すなわち、サポートベクターの集合を$\mathcal{V}$で表現すると、$w$は次で計算できます。
$$
\begin{align}
w&=\sum_{i=1}^{l}\alpha_iy_ix_i\\
&=\sum_{x_i \in \mathcal{V}}\alpha_iy_ix_i
\end{align}
$$
$b$の計算方法
ソフトマージンサポートベクターマシンでも分類超平面の式の$b$の計算方法はハードマージンと同じ方法で求めます。
すなわち、正例、負例のサポートベクター（$0<\alpha_i < C$のとき）をそれぞれ$x_{+}, x_{-}$とし、支持超平面の式に代入すると
$$
\begin{align}
w^Tx_{+}+b&=1,\\
w^Tx_{-}+b&=-1
\end{align}
$$
なので、両式を足し合わせ、$b$に関して解くと
$$
\begin{align}
b&=-\frac{1}{2}w^T(x_{+}+x_{-})\\
&=-\frac{1}{2}\sum_{\mathcal{V}}\alpha_iy_ix_i^T(x_{+}+x_{-})\\
&=-\frac{1}{2}\sum_{\mathcal{V}}(\alpha_iy_ix_i^Tx_{+}+\alpha_iy_ix_i^Tx_{-})
\end{align}
$$
が得られます。
決定関数
サポートベクターのみを使って決定関数$\hat{f}(x)$を表すと
$$
\begin{align}
\hat{f}(x) &= \text{sgn}(w^T x + b)\\
&=\text{sgn}(\sum_{\mathcal{V}}\alpha_iy_ix_i^Tx + b)
\end{align}
$$
これで分類超平面の式が求まりましたので、以上をまとめます。
ソフトマージンサポートベクターマシン
訓練データ
$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}$
が与えられたとき、次の最適化問題を解き、$\alpha$を求めます。
- 目的関数
$$
\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j
$$
を制約条件のもとで$\alpha$に関して最小化
- 制約条件
$$
C \geqq \alpha_i \geqq 0, i=1,\ldots,l,\\
\sum_{i=1}^{l}\alpha_iy_i=0
$$このとき、決定関数は
$$
\hat{f}(x) = \text{sgn}(\sum_{\mathcal{V}}\alpha_iy_ix_i^Tx + b)
$$
ここで、$b$は
$$
\begin{align}
b&=-\frac{1}{2}\sum_{\mathcal{V}}(\alpha_iy_ix_i^Tx_{+}+\alpha_iy_ix_i^Tx_{-})
\end{align}
$$
$\mathcal{V}$はサポートベクターの集合です。
ところで、上のソフトマージンサポートベクターマシンのアルゴリズムを見ると、訓練データ$x_i$は全て内積$x_i^Tx_j$の形でしか現れていないことに気付きます。このことを利用すると次に述べるカーネル法による非線形のサポートベクターマシンが導かれます。
*1:LIBSVM（svm-train）では「-c」オプションを用いてこのパラメータを指定できます。
*2:たとえばLIBSVM（svm-train）のヘルプではサポートベクターマシンの種類を指定する「-s」オプションの説明に「0 -- C-SVC」と書かれています。したがってソフトマージンサポートベクターマシンでデータの分類を学習するのであれば「-s 0」を指定します
サポートベクターマシンとは［最適化問題の解法］
はじめに最適化問題の解法について一般論を述べた後、それをサポートベクターマシンで現れる最適化問題に適用していきます。
最適化問題とは、「ある制約の下で、関数の最小値や最大値を発見すること」で、次のように定式化できます。
なお、目的関数が2次関数で、不等式制約、等式制約は線形関数の最適化問題は2次計画問題(quadratic programming)と呼ばれます。全ての式が線形のときは線形計画問題(linear programming)と呼ばれます。
サポートベクターマシンの目的関数$\frac{||w||^2}{2}$は２次関数、制約条件は線形ですので、サポートベクターマシンの最適化問題は２次計画問題となります。
サポートベクターマシンの最適化問題には等式制約は現れないのですが、ここでは等式制約も含めた最適化問題の一般的な解法であるラグランジュ未定乗数法について述べます。
まず、前述の最適化問題に出てくる目的関数、等式制約、不等式制約を融合した関数$L$（ラグランジュ関数またはラグランジアンと呼ばれます）を定義します。$$
L(z, \alpha, \beta)=f(z)+\sum_{i=1}^{k}\alpha_ig_i(z)+\sum_{i=1}^{m}\beta_ih_i(z)
$$
ここで、
$$
\begin{align}
\alpha&=(\alpha_1,\ldots, \alpha_k)^T, \alpha_i\geqq 0, i=1,\ldots,k\\
\beta&=(\beta_1,\ldots, \beta_k)^T, \beta_i \geqq 0, i=1,\ldots,m
\end{align}
$$
です。この実数ベクトル$\alpha, \beta$はラグランジュ乗数と呼ばれます。
このとき、（制約付き）最適化問題の解は、次のクーン・タッカー(Kuhn-Tucker)の定理を利用して求めることができます。
クーン・タッカー(Kuhn-Tucker)の定理
最適化問題の解が存在する必要十分条件は、その解$x$において次の５つの条件が成り立つような$\alpha, \beta$が存在することです。
$$
\begin{align}
\frac{\partial L(z, \alpha, \beta)}{\partial z}&=0,\\
\frac{\partial L(z, \alpha, \beta)}{\partial \beta}&=0,\\
\alpha_ig_i(z)&=0, i=1,\ldots,k,\\
g_i(z)&\leqq 0, i=1,\ldots,k,\\
\alpha_i&\geqq 0, i=1,\ldots,k,
\end{align}
$$
なお、上の3番目の等式$\alpha_ig_i(z)=0, i=1,\ldots,k$はカルーシュ・クーン・タッカー(Karush-Kuhn-Tucker)の（相補性）条件（KKT条件と略します）として知られています。これから述べるようにKKT条件はサポートベクターマシンの定式化で重要な役割を果たします。
最適化問題を実際に解く前にKKT条件について検討しておきます。
KKT条件から、$\alpha_i$と$g_i(z)$の積は０になるので、次の2つの場合のどちらかが成り立っていることがわかります。
- $g_i(z)< 0$のときは$\alpha_i=0$
- $g_i(z)= 0$のときは$\alpha_i<0$
サポートベクターマシンの不等式制約
$$y(w^Tx+b)\geqq 1$$
にこのKKT条件を適用しましょう。
$z \leftarrow (w, b)$、$g_i(w, b)=1-y_i(w^Tx_i+b)\leqq 0$なのでKKT条件から
- $1-y_i(w^Tx_i+b)< 0$のときは$\alpha_i=0$
- $1-y_i(w^Tx_i+b)= 0$のときは$\alpha_i<0$
のどちらかが成り立ちます。すなわち、
- $y_i(w^Tx_i+b)> 1$のときは$\alpha_i=0$
- $y_i(w^Tx_i+b)= 1$のときは$\alpha_i<0$
となります。
まず、この2番めのケースに注目すると、この式ははまさしく$x_i$が支持超平面（$y(w^Tx+b)=1$）にちょうど乗っている、すなわち$x_i$がサポートベクターである場合に相当します。したがって、$x_i$がサポートベクターであるときは、対応するラグランジュ乗数$\alpha_i<0$となります。
つぎに、1番目のケースを見ます。$x_i$が支持超平面の正例側または負例側のどちらかに存在する、すなわち$x_i$はサポートベクターではないときは、ラグランジュ乗数$\alpha_i=0$となることを意味しています。
以上をまとめると、最適化問題の解のうち、$\alpha_i<0$に対応するデータ$x_i$はサポートベクターとなります。
さらに、サポートベクターではないデータ$x_i$では$\alpha_i=0$となることを利用すれば分類超平面の$w$の計算回数を削減できます（後述）。
それではクーン・タッカー(Kuhn-Tucker)の定理を使って最適化問題を実際に解いていきます。
まずサポートベクターマシンの最適化問題に対するラグランジュ関数$L(w, b, \alpha, \beta)$を求めます。
$$
\begin{align}
f(w, b)&=\frac{||w||^2}{2},\\
g_i(w, b)&=1-y_i(w^Tx_i+b), i=1,\ldots,l,\\
h_i(w, b)&=0, i=1,\ldots,m
\end{align}
$$
とおけるので、このときラグランジュ関数は、
$$
\begin{align}
L(w, b, \alpha, \beta) &= f(w, b)+\sum_{i=1}^{l}\alpha_ig_i(w, b)+\sum_{i=1}^{m}\beta_ih_i(w, b)\\
&= \frac{||w||^2}{2}+\sum_{i=1}^{l}\alpha_i(1-y_i(w^Tx_i+b))
\end{align}
$$
となり、これを$w$、$b$、$\beta$で偏微分すると
$$
\begin{align}
\frac{\partial L(w, b, \alpha, \beta)}{\partial w}&=w-\sum_{i=1}^{l}\alpha_iy_ix_i\\
&=0,\\
\frac{\partial L(w, b, \alpha, \beta)}{\partial b}&=-\sum_{i=1}^{l}\alpha_iy_i\\
&=0\\
\frac{\partial L(w, b, \alpha, \beta)}{\partial \beta}&=0
\end{align}
$$
ここで$\partial ||w||^2/\partial w = \partial w^Tw/\partial w =2w$であることを使いました。
整理すると、
$$
\begin{align}
w&=\sum_{i=1}^{l}\alpha_iy_ix_i,\\
\sum_{i=1}^{l}\alpha_iy_i&=0
\end{align}
$$
となります。
しかしこれらの式だけでは$\alpha$を決めることができません。
そこでこれまでの最適化問題とは別の$\alpha$に関する最適化問題（双対問題）を導き、これを解くことにより$\alpha$を求めます。
$w=\sum_{i=1}^{l}\alpha_iy_ix_i$なので、これをラグランジュ関数$L(w, b, \alpha, \beta)$に代入した関数を
$$
L_D(\alpha, \beta)=L(\sum_{i=1}^{l}\alpha_iy_ix_i, b, \alpha, \beta)
$$
と定義します。
最適化問題（双対問題）
- 目的関数
$$L_D(\alpha, \beta)$$
を制約条件のもとで$\alpha$に関して最小化
- 制約条件
$$\alpha_i\geqq 0, i=1,\ldots,l\\
\sum_{i=1}^{l}\alpha_iy_i=0$$
この最適化問題は主問題に対する双対問題と呼ばれます。なお、制約条件に$\beta$は現れないことに注意してください。
ラグランジュ関数$L_D(b, \alpha, \beta)$を具体的に計算しましょう。
$$
\begin{align}
L_D(\alpha, \beta)&=L(\sum_{i=1}^{l}\alpha_iy_ix_i, b, \alpha, \beta)\\
&=\frac{||w||^2}{2}+\sum_{i=1}^{l}\alpha_i(1-y_i(w^Tx_i+b))\\
&=\frac{w^Tw}{2}+\sum_{i=1}^{l}\alpha_i(1-y_i(w^Tx_i+b))\\
&=\frac{(\sum_{i=1}^{l}\alpha_iy_ix_i)^T(\sum_{i=1}^{l}\alpha_iy_ix_i)}{2}+\sum_{i=1}^{l}\alpha_i\{1-y_i((\sum_{i=1}^{l}\alpha_iy_ix_i)^Tx_i+b)\} \\
&=\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j+\sum_{i=1}^{l}\alpha_i-\sum_{i=1}^{l}\alpha_iy_i(\sum_{i=1}^{l}\alpha_iy_ix_i)^Tx_i-b(\sum_{i=1}^{k}\alpha_iy_i)\\
&=\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j+\sum_{i=1}^{l}\alpha_i-\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j (\because \sum_{i=1}^{l}\alpha_iy_i=0)\\
&=\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j
\end{align}
$$
したがって、双対問題の目的関数：
$$
\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j
$$
を制約条件：
$$\alpha_i\geqq 0, i=1,\ldots,l\\
\sum_{i=1}^{l}\alpha_iy_i=0$$
のもとで最小化することにより$\alpha$が求めれば良いことになります。
$\alpha$がわかれば、$w$も$w=\sum_{i=1}^{l}\alpha_iy_ix_i$から求まります。
この時点で分類超平面の式のパラメータ$w$は求まりましたが、まだ$b$がわかりません。
ここでサポートベクターが支持超平面（$y(w^Tx+b)=1$）に乗っているという事実を思い出せば、この式から$b$が直接計算できることが分かります。すなわち、$x_{+}, x_{-}$をそれぞれ正例、負例のサポートベクターとすれば、これを支持超平面の式に代入して
$$
\begin{align}
w^Tx_{+}+b&=1,\\
w^Tx_{-}+b&=-1
\end{align}
$$
が得られ、この２式を足し合わせて$b$に関して解くと、$b$を求める式
$$
b=-\frac{1}{2}w^T(x_{+}+x_{-})
$$
が得られます。
ところで以前、KKT条件からサポートベクターではないデータでは$\alpha_i=0$であることを見ました。サポートベクターとなるデータは訓練データ全体と比べると圧倒的に少ないので、これを利用して分類超平面の$w$の計算量を削減できます。
すなわち、サポートベクターの集合を$\mathcal{V}$で表現すると、$w$は次式で計算できます。
$$
\begin{align}
w&=\sum_{i=1}^{l}\alpha_iy_ix_i\\
&=\sum_{x_i \in \mathcal{V}}\alpha_iy_ix_i
\end{align}
$$
これで分類超平面を求める式が出揃いました。以上をまとめるとつぎのようになります。
ハードマージンサポートベクターマシン
線形分離可能な訓練データ
$$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times {-1, 1}$$
が与えられたとき、次の最適化問題を解き$\alpha$を求めます。
- 目的関数
$$
\sum_{i=1}^{l}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{l}\alpha_i\alpha_jy_iy_jx_i^Tx_j
$$
を次の最適条件のもとで$\alpha$に関して最小化
- 制約条件
$$\alpha_i\geqq 0, i=1,\ldots,l$$
$$\sum_{i=1}^{l}\alpha_iy_i=0$$このとき決定関数$\hat{f}(x)$は
$$
\hat{f}(x) = \text{sgn}(w^T x + b)
$$
です。ここで、$w,b$は
$$
w=\sum_{\mathcal{V}}\alpha_iy_ix_i\\
b=-\frac{1}{2}w^T(x_{+}+x_{-})
$$
で与えられます（$x_{+}, x_{-}$はそれぞれ正例、負例のサポートベクター）。
これで、サポートベクターマシンの基礎が一応完成しましました。
とくに訓練データが線形分離可能な場合のアルゴリズムのことをハードマージンサポートベクターマシンと呼ぶことがあります。
残念ながら、ハードマージンサポートベクターマシンを現実の問題に適用しても高い分類性能が得られない可能性が高いと考えられます。なぜなら、現実の状況では正確なデータだけでなく測定誤差や人為的ミスなどによるデータも発生するので、当然訓練データにもそのようなデータが含まれています。したがって、訓練データが線形分離可能であることはあまり期待できないため、ハードマージンサポートベクターマシンの前提が満たされません。
また、支持超平面から離れたデータが誤差を含んでいてもこれらのデータが支持超平面の決定に影響を与えることはないので問題はありませんが、クラス境界近くにある誤差データは分類超平面の決定に直接的に影響するため精度高い分類性能は期待できなさそうです。
このように誤差を含むデータの影響を少なくし、サポートベクターマシン本来の高い分類性能を発揮できるように、サポートベクターマシンにスラック変数を導入します。これが次に説明するソフトマージンサポートベクターマシンです。
サポートベクターマシンとは［ハードマージンサポートベクターマシン］
まずはじめに訓練データが線形分離可能な場合について定式化します。この場合のサポートベクターマシンをハードマージンサポートベクターマシンと呼びます。
線形分離できないデータへの拡張（ソフトマージンサポートベクターマシン）については次の記事で説明します。
準備
まず初めに記号を定義しておきます。
訓練データのサイズ（個数）：$l$
訓練データの次元：$n$
とすると、属性データの集合$X$とそれに対応するクラスの集合$Y$をそれぞれ次のように表記します。
$$
\begin{align}
X &= \{x_1, x_2, \ldots, x_l\}, x_i \in \Re^n\\
Y &= \{y_1, y_2, \ldots, y_l\}, y_i \in \{-1, 1\}
\end{align}
$$
$X$の各要素$x_i$が１個の訓練データに相当し、それぞれ$n$次元ベクトルであることに注意してください。
またデータ$x_i$のクラスは$y_i$となります。後の定式化が楽になるようにクラスの値は１（正例）または－１（負例）である仮定します。
たとえば、スパム判定問題の場合には、メール$x_i$がスパムであるときそのクラス$y_i$は１，スパムでないときは－１となります。
このとき訓練データの集合$D$は
$$
\begin{align}
D &= \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}\\
&= X \times Y
\end{align}
$$
と表現できます。
分類超平面と決定関数
さらに両クラスのデータを２つに綺麗に分類（分割）する超平面（分類超平面と呼びます）の方程式を
$$
w^Tx+b=0
$$
と表記します。ここで$w \in \Re^n$は超平面から直角に正例データのある方向に向かうベクトル（法線ベクトル）、$b \in \Re$は切片に相当します。また、$w^Tx$は$w$と$x$の内積、$x^T$の$T$は転置記号です。
念の為に書いておきますと、この式の$x$は単に$\Re^n$空間の変数を表しているだけです。訓練データの$x_i$ではありません。
この分類超平面の一方に正例データ、反対側に負例データがありますので、分類超平面の式を求めさえすれば、$x$に未知のデータを代入した結果をもとにこのデータを正例または負例に分類することができます。
すなわち、サポートベクターマシンで分類超平面の式の$w$と$b$が求まれば、未知データ$x$の属するクラス$y$は次の式（決定関数、識別関数と呼ばれます）で求めることができます。
$$
\begin{align}
y &=
\begin{cases}
1, & \text{if } w^T x +b \geqq 0 \\
-1, & \text{if } w^T x +b < 0
\end{cases}\\
&= \text{sgn}(w^T x + b)
\end{align}
$$
ここで$\text{sgn}(x)$は符号関数で次で定義されます。
$$
\text{sgn}(x) =
\begin{cases}
1, & \text{if } x \geqq 0 \\
-1, & \text{if } x < 0
\end{cases}
$$
このように未知のデータでも精度良く分類できる分類超平面（式の$w$と$b$）を求めることがサポートベクターマシンアルゴリズムの目的となります。
これからサポートベクターマシンを構成していくわけですが、その前に用語を定義します。
- 支持超平面
- 分類超平面から最短距離にある正例または負例のデータを通り、かつ、分類超平面と平行な超平面のこと。下図の点線が支持超平面となります。
- マージン
- 分類超平面を挟む２枚の支持超平面の間の距離のこと（下図参照）。
- サポートベクター
- 支持超平面の上にちょうど乗っているデータのこと（下図参照）。
さて、SVMでは分類超平面を求めるわけですが、いまのままでは分類超平面が一つに決まりません。実際、正例と負例を分ける超平面は無数に存在します。たとえば、上の図にあるように緑の実線で表された超平面も、黒の実線の超平面も共に分類超平面の候補となります。
このままでは困るので、サポートベクターマシンでは、分類性能の高い分類超平面を一意に決めるために、超平面に次の条件を課します。
- 分類超平面は正例と負例のちょうど真ん中に位置する
- マージンを最大にする
したがってこれらの条件を満たす分類超平面を求めることになるわけですが、上図を見ればわかるように、実際は少数のサポートベクターの情報さえあれば分類超平面を決定できます。それ以外の大部分の訓練データは分類超平面の決定に影響しないということになります。
分類超平面の導出
支持超平面は分類超平面（$w^Tx+b=0$）に平行なので、正例側および負例側の支持超平面は、それぞれ次のように表現できます（$k>0$）。
$$
\begin{align}
w^Tx+b &= k\\
w^Tx+b &= -k
\end{align}
$$
直線（平面）の方程式を定数倍しても、元の直線（平面）と同じなので、支持超平面の式（$w^Tx+b=\pm k$）を$1/k$倍して
$$\frac{1}{k}w^Tx+\frac{1}{k}b=\pm 1$$
とできます。$1/k$があると式変形が煩雑になるため、
$$\frac{1}{k}w \rightarrow w, \frac{1}{k}b \rightarrow b$$
と書き直すことにより、最終的に支持超平面は
$$w^Tx+b=\pm 1$$
と表現できます。同様に分類超平面はつぎの式になります。
$$w^Tx+b=0$$
では、分類超平面の$w$と$b$を具体的に求めましょう。
まずは、上述した分類超平面の条件である、マージンの最大化、を検討します。
$x_p$、$x_q$をそれぞれ正例側、負例側のサポートベクターとすると、支持超平面間のマージン$m$はつぎのようになります（上図を参照）。
$$
\begin{align}
m &= ||x_p-x_q||\cos\theta \\
&=||x_p-x_q||\frac{w^T(x_p-x_q)}{||w|| ||x_p-x_q||} (\because a^Tb=||a|| ||b||\cos\theta)\\
&=\frac{w^Tx_p-w^Tx_q}{||w||}\\
&=\frac{(1-b)-(-1-b)}{||w||} (\because w^Tx_p+b=1, w^Tx_q+b=-1)\\
&=\frac{2}{||w||}
\end{align}
$$
マージン$m$の最大化は、
$$
\begin{align}
\max_{w, b}\frac{2}{||w||} &\Rightarrow \min_{w, b} 2||w||\\
&\Rightarrow \min_{w, b} 2||w||^2\\
&\Rightarrow \min_{w, b} \frac{||w||^2}{2}
\end{align}
$$
と表現出来ます。ここで、$\max_{w, b}, \min_{w, b}$はそれぞれ$w$と$b$に関する最大化、最小化を意味する記号です。また最小化（最大化）の対象となる式（$||w||^2/2$）のことを目的関数と呼びます。
上の式変形についてそれぞれコメントします。
- $1/||w||$を最大化するということは、$||w||$を最小化することと同じ
- $||w||$を最小化するということは、その2乗$||w||^2$を最小化することと同じ
- $||w||^2$を最小化することは、$||w||^2/2$を最小化することと同じ（後の式変形の都合上２で割っています）
以上より、マージンを最大化するには、目的関数$||w||^2/2$を最小化する$w$を求めれば良いことが分かりました。
この目的関数$||w||^2/2$の式だけをみると$w=0$とすれば最小化できますが、実際には後述する制約条件（正例・負例は分類超平面の片側のみに存在という条件）を考慮に入れて最小化しなければいけないため、$w=0$とすることはできません。
つぎに、もう一つの分類超平面の条件である、分類超平面は正例と負例のちょうど中間に位置、について考えます。後述しますが、この条件から導かれる関係式が上述の最小化問題の制約条件になります。
分類超平面は正例と負例のちょうど中間に位置することから、分類超平面の片側には正例データまたは負例データのみ存在します。これを数式で表現することを考えます（下の図を参照）。
まず正例側にある支持超平面の場合から。
図の$x_A, c, z$はすべてベクトルなので、$x_A=c+z$が成り立ちます。この式と$w$の内積を計算すると
$$
\begin{align}
w^Tz&=w^Tx_A-w^Tc\\
&=w^Tx_A-(1-b) (\because w^Tc+b=1)
\end{align}
$$
正例データは支持超平面（$w^Tc+b=1$）の左側にのみある、すなわち角度$\theta$は$-\frac{\pi}{2} \leqq \theta \leqq \frac{\pi}{2}$なので
$$
\begin{align}
w^Tz&=||w|| ||z||\cos\theta\\
&\geqq 0
\end{align}
$$
となります。これらの式から
$$
\begin{align}
w^Tx_A-(1-b)\geqq 0\\
\Downarrow\\
w^Tx_A+b\geqq 1
\end{align}
$$
が成り立ちます。
負例データの場合も同様の方法で求めます。$x_B=c+z$と$w$の内積を計算すると
$$
\begin{align}
w^Tz&=w^Tx_B-w^Tc\\
&=w^Tx_B-(-1-b) (\because w^Tc+b=-1)
\end{align}
$$
負例データは支持超平面（$w^Tc+b=-1$）の右側にのみ存在、すなわち角度$\theta$は$\frac{\pi}{2} \leqq \theta \leqq \frac{3\pi}{2}$なので
$$
\begin{align}
w^Tz&=||w|| ||z||\cos\theta\\
& \leqq 0
\end{align}
$$
となります。これらの式から
$$
\begin{align}
w^Tx_B-(-1-b)\leqq 0\\
\Downarrow\\
w^Tx_B+b\leqq -1
\end{align}
$$
が成り立ちます。
以上をまとめておきます。
正例データのクラスは１（$y=1$）、負例データのクラスは－１（$y=-1$）であったことを思い出すと、あるデータ$x$のクラスは
$$
y =
\begin{cases}
1, & \text{if } w^T x +b \geqq 1 \\
-1, & \text{if } w^T x +b < -1
\end{cases}
$$
で求められます。
この式のままですと場合分けの必要があるため取り扱いが面倒なので一つの式で表現することを考えます。
正例データは$y=1$なので、上の正例の場合の式を$y$を用いて書き換えると
$$
w^Tx+b\geqq 1\\
\Downarrow \\
1(w^Tx+b)\geqq 1\\
\Downarrow\\
y(w^Tx+b)\geqq 1
$$
となります。同様に負例データのクラスは$y=-1$なので
$$
w^Tx+b\leqq -1\\
\Downarrow \\
-1(w^Tx+b)\geqq 1\\
\Downarrow \\
y(w^Tx+b)\geqq 1
$$
となります。これで正例・負例の場合をまとめて一つの式
$$
y(w^Tx+b)\geqq 1
$$
で表現できました。
これが前述した目的関数$||w||^2/2$を最小化する際の制約条件となります。
これで必要な関係式が揃いましたので、分類超平面を求める問題は次の最適化問題として定式化できます。
最大マージン分類器
線形分離可能な訓練データ
$$D = \{(x_1, y_1), \ldots, (x_l, y_l)\}, (x_i, y_i) \in \Re^n \times \{-1, 1\}$$
が与えられたとき、次の制約付き最適化問題を解き最大マージン分類超平面$w^Tx+b=0$を求めます。
- 目的関数（コスト関数）の最小化
$$\min_{w, b} \frac{||w||^2}{2}$$
- 制約条件
$$y(w^Tx+b)\geqq 1$$
これからこの最適化問題を解いていきます。
サポートベクターマシンとは［はじめに］
本記事ではサポートベクターマシンについて説明します。Wikipediaによるとサポートベクターマシン（support vector machine、SVM）とは
教師あり学習を用いる識別手法の一つである。パターン認識や回帰分析へ適用できる。
サポートベクターマシンは、現在知られている多くの手法の中で一番認識性能が優れた学習モデルの一つである。サポートベクターマシンがすぐれた認識性能を発揮することができる理由は、未学習データに対して高い識別性能を得るための工夫があるためである。
Wikipediaで言及されている、未知データに対する識別能力のことを汎化能力といいます。実際、数ある識別アルゴリズムのなかでも、サポートベクターマシンの汎化能力はトップクラスです。
サポートベクターマシンは分類問題・回帰問題・新規性発見などの問題に適用できますが、問題の種類ごとに定式化の方法が若干異なります。
この記事ではこのなかでもっとも基本的な分類問題に対するサポートベクターマシンについて述べることにします。これで原理が理解できれば、回帰問題や新規性発見でのアルゴリズム理解はそれほど難しくないと思います。
さて、サポートベクターマシンのアルゴリズムの流れを非常に単純に表現すると次のようになります。
- 学習フェーズ：事前に準備された訓練データ（事前に正解がわかっているデータ）を使って、高い分類性能が得られるようにサポートベクターマシンに学習させます
- 予測フェーズ：未知データを学習済みサポートベクターマシンに与えると、そのデータがどこに分類されるのか、教えてくれます
これらを詳細について説明する前に、ここでサポートベクターマシンの特徴を簡単にまとめておきます。
- ２値分類学習器
サポートベクターマシンは、分類対象が２種類だけの２値データ（２クラスデータ、２群データともいいます）に対する分類器です。したがって、２クラスに分類されるデータを与えると、そのデータをどちらかのクラスに分類します。
たとえば、スパム判定問題では、メールがスパムなのか、そうでないかの２種類しかありませんので、サポートベクターマシンを適用してメールのスパム判定が可能になります。
ところで、２クラスからなる分類対象うち、注目しているクラスに属すデータを正例、それ以外のデータを負例を呼ぶことがあります。スパム判定問題の場合でいうと、正例はスパムメールデータ、負例はスパムでないメールデータ、となります。
実は代表的なサポートベクターマシンソフトウェアを使うと、３クラス以上の多クラスデータも分類できます。上述したようにサポートベクターマシンは２クラスデータしか分類できませんので、多クラスデータに適用するには工夫が必要です。実際ソフトウェア内部では、２クラスデータ用のサポートベクターマシンを多数生成することにより多クラスデータを分類しています。ペアワイズ法やone-versus-rest法が用いられますが、詳細は参考文献を参照してください。
- 線形分離可能なデータに対する分類器
線形分離可能とは、$n$次元データが与えられたとき、正例のデータと負例のデータを$n-1$次元超平面で綺麗に２つに分けられることを言います。
たとえば、２次元の２値データを直線で綺麗に正例・負例の２つに分離できればそのデータは線形分離ということになります。サポートベクターマシンは線形分離可能なデータに対する分類器です。この分類超平面は無数に存在しますが、サポートベクターマシンでは正例・負例データとの距離が最大になるように分類超平面を選びます。
それでは２値クラス分類問題に対するサポートベクターマシンを定式化していきます。
Rの基本グラフィックス機能またはggplot2を使って地図を描くには
Rに元から備わっているグラフィックス機能とその機能を拡張するggplot2で日本地図を表示する方法について説明します。
地図データの準備
日本地図のシェープファイルを入手します。入手先は２ヶ所（Global Administrative AreasおよびESRIジャパン株式会社）ありますので、順に説明します。なお、両者のデータの中身は異なるようなので注意してください。
Global Administrative Areasからダウンロードするには、Downloadにアクセスして、プルダウンメニューの「Country」で「Japan」を選択（「File format」メニューは「Shapefile」のままで）すると日本地図のシェープファイルのダウンロードページにジャンプします。
ちなみに、「File format」メニューで「R (SpatialPolygonsDataFrame)」を選択すると、バイナリ書式ファイルデータ（.Rdata）をダウンロードすることができます。好みに応じで選択してください。
ESRIジャパン株式会社からダウンロードするには、全国市区町村界データにアクセス後、指示に従ってダウンロードして下さい。
ところで、全国市区町村界データのReadme.txtには使用規定が書かれており、そこに次の条項あります。
第2条 遵守事項
1. お客様は、本データをArcGISシリーズ（アメリカ合衆国カリフォルニア州法人である
Environmental Systems Research Institute, Inc.の地理情報システムに関する
一連のソフトウェア製品）以外のソフトウェアで使用することはできません。
これを読むとRでは使用できないと思われるのですが、RjpWiki：ShapeFileライブラリの「ESRIジャパン社が公開している全国市区町村界データの加工・再公開の許可について」を読むとRでも使用できそうにも思えるのですが、どうなのでしょうか。
というわけで、以降ではGlobal Administrative Areasのシェープファイルを用いて日本地図を表示します。
データをダウンロードしたら適当なディレクトリにファイルを解凍してください。
Rライブラリの準備
地図を表示するには次のパッケージが必要となりますので、事前にインストールしておきます。
シェープファイルの読み込み
Global Administrative Areasのデータを解凍すると分かりますが、次の３つの異なるレベルのシェープファイルが含まれています。
- JPN_adm0.shp：日本列島の全体地図
- JPN_adm1.shp：都道府県地図
- JPN_adm2.shp：市区町村地図
たとえば、JPN_adm1.shp（都道府県地図）を読み込むには次のようにします。
library(maptools) library(spsurvey) jpn <- readShapePoly("JPN_adm1.shp") # シェープファイルはカレントディレクトにあるとします jpn.dbf <- read.dbf("JPN_adm1.dbf") # dbfファイルの読み込み
readShapePolyはmaptoolsパッケージの関数です。読み込み時間は１秒くらいでした。
spsurveyパッケージのread.shape関数を使ってもシェープファイルを読み込めますが、非常に時間がかかりますので、あまりオススメはできません。実際、都道府県シェープファイルの読み込みに６５秒、市区町村で７３２秒かかりました（もちろん環境依存）。
dbfファイルの中身を表示してみます。下記の情報を使うと一部の地域のみ表示することも可能になります（例は後述します）。
jpn.dbf ID_0 ISO NAME_0 ID_1 NAME_1 VARNAME_1 NL_NAME_1 HASC_1 CC_1 TYPE_1 ENGTYPE_1 VALIDFR_1 VALIDTO_1 REMARKS_1 Shape_Leng Shape_Area area_mdm 1 116 JPN Japan 1457 Aichi Aiti ??? JP.AI <NA> Ken Prefecture Unknown Unknown <NA> 8.194165 0.5090549 0.5090549 2 116 JPN Japan 1458 Akita <NA> ??? JP.AK <NA> Ken Prefecture Unknown Unknown <NA> 8.995578 1.2236946 1.2236946 3 116 JPN Japan 1459 Aomori <NA> ??? JP.AO <NA> Ken Prefecture Unknown Unknown <NA> 10.431820 1.0315494 1.0315494 4 116 JPN Japan 1460 Chiba Tiba|Tsiba ??? JP.CH <NA> Ken Prefecture Unknown Unknown <NA> 6.571243 0.5011059 0.5011059 5 116 JPN Japan 1461 Ehime <NA> ??? JP.EH <NA> Ken Prefecture Unknown Unknown <NA> 15.999670 0.5586251 0.5586251 6 116 JPN Japan 1462 Fukui Hukui ??? JP.FI <NA> Ken Prefecture Unknown Unknown <NA> 6.790149 0.4224879 0.4224879 7 116 JPN Japan 1463 Fukuoka Hukuoka ??? JP.FO <NA> Ken Prefecture Unknown Unknown <NA> 8.266049 0.4862017 0.4862017 8 116 JPN Japan 1464 Fukushima Hukusima ??? JP.FS <NA> Ken Prefecture Unknown Unknown <NA> 7.533269 1.3975740 1.3975740 9 116 JPN Japan 1465 Gifu Gihu ??? JP.GF <NA> Ken Prefecture Unknown Unknown <NA> 6.635546 1.0570861 1.0570861 10 116 JPN Japan 1466 Gunma GunmaGumma ??? JP.GM <NA> Ken Prefecture Unknown Unknown <NA> 4.920751 0.6406968 0.6406968 11 116 JPN Japan 1467 Hiroshima Hirosima ??? JP.HS <NA> Ken Prefecture Unknown Unknown <NA> 13.029234 0.8284376 0.8284376 12 116 JPN Japan 1468 Hokkaido Ezo|Yeso|Yezo ??? JP.HK <NA> Do Circuit Unknown Unknown <NA> 35.144123 8.6717107 8.6717107 13 116 JPN Japan 1469 Hyogo Hiogo ??? JP.HG <NA> Ken Prefecture Unknown Unknown <NA> 9.991132 0.8280074 0.8280074 14 116 JPN Japan 1470 Ibaraki <NA> ??? JP.IB <NA> Ken Prefecture Unknown Unknown <NA> 5.595910 0.6125375 0.6125375 15 116 JPN Japan 1471 Ishikawa Isikawa ??? JP.IS <NA> Ken Prefecture Unknown Unknown <NA> 7.145355 0.4275550 0.4275550 16 116 JPN Japan 1472 Iwate <NA> ??? JP.IW <NA> Ken Prefecture Unknown Unknown <NA> 9.800987 1.5975500 1.5975500 17 116 JPN Japan 1473 Kagawa <NA> ??? JP.KG <NA> Ken Prefecture Unknown Unknown <NA> 5.257244 0.1689363 0.1689363 18 116 JPN Japan 1474 Kagoshima <NA> ???? JP.KS <NA> Ken Prefecture Unknown Unknown <NA> 25.913128 0.8627851 0.8627851 19 116 JPN Japan 1475 Kanagawa <NA> ???? JP.KN <NA> Ken Prefecture Unknown Unknown <NA> 5.303150 0.2365421 0.2365421 20 116 JPN Japan 1476 Kochi Koti ??? JP.KC <NA> Ken Prefecture Unknown Unknown <NA> 9.432022 0.6827940 0.6827940 21 116 JPN Japan 1477 Kumamoto <NA> ??? JP.KM <NA> Ken Prefecture Unknown Unknown <NA> 13.608921 0.7266430 0.7266430 22 116 JPN Japan 1478 Kyoto Kioto ??? JP.KY <NA> Fu Urban Prefecture Unknown Unknown <NA> 6.139137 0.4578284 0.4578284 23 116 JPN Japan 1479 Mie Miye ??? JP.ME <NA> Ken Prefecture Unknown Unknown <NA> 12.134027 0.5670699 0.5670699 24 116 JPN Japan 1480 Miyagi <NA> ??? JP.MG <NA> Ken Prefecture Unknown Unknown <NA> 11.110625 0.7453862 0.7453862 25 116 JPN Japan 1481 Miyazaki <NA> ??? JP.MZ <NA> Ken Prefecture Unknown Unknown <NA> 7.626958 0.7343826 0.7343826 26 116 JPN Japan 1482 Nagano <NA> ??? JP.NN <NA> Ken Prefecture Unknown Unknown <NA> 7.662159 1.3596386 1.3596386 27 116 JPN Japan 1483 Naoasaki Nagasaki ??? JP.NS <NA> Ken Prefecture Unknown Unknown <NA> 35.052216 0.4067374 0.4067374 28 116 JPN Japan 1484 Nara <NA> ??? JP.NR <NA> Ken Prefecture Unknown Unknown <NA> 3.485988 0.3607166 0.3607166 29 116 JPN Japan 1485 Niigata <NA> ??? JP.NI <NA> Ken Prefecture Unknown Unknown <NA> 11.208076 1.2932382 1.2932382 30 116 JPN Japan 1486 Oita <NA> ??? JP.OT <NA> Ken Prefecture Unknown Unknown <NA> 9.383343 0.6134952 0.6134952 31 116 JPN Japan 1487 Okayama <NA> ??? JP.OY <NA> Ken Prefecture Unknown Unknown <NA> 10.899804 0.7159576 0.7159576 32 116 JPN Japan 1488 Okinawa <NA> ??? JP.ON <NA> Ken Prefecture 19720515 Unknown <NA> 17.981032 0.2105791 0.2105791 33 116 JPN Japan 1489 Osaka <NA> ??? JP.OS <NA> Fu Urban Prefecture Unknown Unknown <NA> 5.167053 0.1880098 0.1880098 34 116 JPN Japan 1490 Saga <NA> ??? JP.SG <NA> Ken Prefecture Unknown Unknown <NA> 4.872416 0.2376768 0.2376768 35 116 JPN Japan 1491 Saitama <NA> ??? JP.ST <NA> Ken Prefecture Unknown Unknown <NA> 3.194298 0.3793111 0.3793111 36 116 JPN Japan 1492 Shiga Siga ??? JP.SH <NA> Ken Prefecture Unknown Unknown <NA> 3.300463 0.3976926 0.3976926 37 116 JPN Japan 1493 Shimane Simane ??? JP.SM <NA> Ken Prefecture Unknown Unknown <NA> 11.457792 0.6690916 0.6690916 38 116 JPN Japan 1494 Shizuoka Sizuoka ??? JP.SZ <NA> Ken Prefecture Unknown Unknown <NA> 13.495757 0.7898159 0.7898159 39 116 JPN Japan 1495 Tochigi Totigi ??? JP.TC <NA> Ken Prefecture Unknown Unknown <NA> 3.694183 0.6456594 0.6456594 40 116 JPN Japan 1496 Tokushima Tokusima ??? JP.TS <NA> Ken Prefecture Unknown Unknown <NA> 5.862819 0.4024895 0.4024895 41 116 JPN Japan 1497 Tokyo Edo|Yedo|Tokio|T｢quio ??? JP.TK <NA> To Metropolis Unknown Unknown <NA> 4.466625 0.1793922 0.1793922 42 116 JPN Japan 1498 Tottori <NA> ??? JP.TT <NA> Ken Prefecture Unknown Unknown <NA> 4.444174 0.3450015 0.3450015 43 116 JPN Japan 1499 Toyama <NA> ??? JP.TY <NA> Ken Prefecture Unknown Unknown <NA> 4.233530 0.4302356 0.4302356 44 116 JPN Japan 1500 Wakayama <NA> ???? JP.WK <NA> Ken Prefecture Unknown Unknown <NA> 7.073690 0.4533757 0.4533757 45 116 JPN Japan 1501 Yamagata <NA> ??? JP.YT <NA> Ken Prefecture Unknown Unknown <NA> 5.430651 0.9657578 0.9657578 46 116 JPN Japan 1502 Yamaguchi Yamaguti ??? JP.YC <NA> Ken Prefecture Unknown Unknown <NA> 15.226488 0.5984350 0.5984350 47 116 JPN Japan 1503 Yamanashi Yamanasi ??? JP.YN <NA> Ken Prefecture Unknown Unknown <NA> 3.404750 0.4447103 0.4447103
ちなみに、バイナリ書式ファイルデータ（.Rdata）をダウンロードした場合は、load関数でデータを読み込みます。
load("JPN_adm1.RData") # オブジェクトgadmにシェープファイルデータが格納されています。
地図の表示
基本グラフィックスまたはggplot2を使って地図を表示します。
地図の表示範囲を指定できると便利なので、必要に応じて次のように表示範囲の経度・緯度を設定します。
- 日本列島全体
xlim <- c(122, 154) # 経度 ylim <- c(20, 46) # 緯度
この設定ですと日本の領土全体を含むので表示範囲が広すぎるかもしれません。
- 日本列島（沖縄を含める）
xlim <- c(126, 146) # 経度 ylim <- c(26, 46) # 緯度
- 日本列島（沖縄を含めない）
xlim <- c(128, 146) # 経度 ylim <- c(30, 46) # 緯度
都道府県地図を表示する場合を想定して、県ごとの色を適当に決めてから基本グラフィックスで地図を表示しましょう。
library(RColorBrewer) col <- sample(1:8,size=47,replace=TRUE) # 県ごとの色 plot(jpn, xlim=xlim, ylim=ylim, col=brewer.pal(8,"Accent")[col]) # 地図をプロット。色はAccentパレットを使用。
すると次のような都道府県地図と市区町村地図（JPN_adm2.shpを読み込んだ場合）が表示されます。
もし、東京都だけを表示したければ、シェープファイルデータを格納している変数（ここの場合jpn）から東京都のデータだけを抜き出してプロットすれば欲しい地図が得られます。
上の方で表示したdbfファイルのNAME_1カラムに県名が入っていますので、それを利用して抽出します。
東京都に対応するNAME_1は「Tokyo」です。ちなみに長崎県は「Naoasaki」（Nagasakiではない）となっているので注意してください。
次のスクリプトは東京都の地図に新宿区役所の場所を表示するためのサンプルです。
plot(jpn[jpn$NAME_1=="Tokyo",]) # 東京都のデータだけを抽出。JPN_adm2.shpを読み込んでいるとします。 points(139.7036, 35.69389, lwd=20, col="red") # 新宿区役所の位置に打点 text(139.7036, 35.69389, "新宿", col="blue", adj = c(-0.3,0.5), cex=2) # 文字を表示
次にggplot2を使って地図を表示しましょう。
ggplot2で地図を表示するときの注意点として、事前にfortify関数でシェープファイルデータをggplot2に合致するデータ形式に変換しておく必要があります。
Rスクリプト例と結果を示します。
library(ggplot2) library(maptools) jpn <- readShapePoly("JPN_adm1.shp") # シェープファイルはカレントディレクトにあるとします map <- fortify(jpn) # 日本全体 #map <- fortify(jpn[jpn$NAME_1=="Tokyo",]) # 東京都だけ表示したい場合 xlim <- c(128, 146) # 経度 ylim <- c(30, 46) # 緯度 col <- sample(1:8,size=47,replace=TRUE) # 県ごとの色 ggplot() + geom_polygon(aes(long,lat,group=group,fill=as.character(col[as.integer(id)+1])),color="black",data=map) + coord_fixed(ratio=1) + geom_point(aes(x=139.7036, y=35.69389),color="blue",size=10) + geom_text(aes(x=139.7036, y=35.69389),label="新宿",hjust=-0.2,color="red",size=10) + guides(col=FALSE, fill=FALSE) + xlim(xlim) + ylim(ylim) + scale_fill_brewer(type="qual",palette="Accent")
参考サイト：
「金塊か、キノコ料理か」（外れ値検出問題）を解く［おわりに］
さまざまな外れ値検出法を用いて、100個の玉の中から貴重な石を一つだけ選び出す課題を解いてきました。
手法により結果は異なりますが、８７番目のデータが外れ値である可能性が最も高そうです。
最後にこの出題自体に対するコメントを述べたいと思います。
まず、外れ値の意味ですが、Wikipediaによると
統計において他の値から大きく外れた値
のことですし、これとは別の定義として、次のホーキンス（Hawkins）による定義があります。
ほかの測定値から大きく外れているため、異なるメカニズムで生成されたとの疑いを抱かせるような観測値のこと
このように、外れ値は定性的に定義されるだけで、定量的には定義されていませんし、定量的かつ一般的に定義すること自体が困難です。
このような事情のため、外れ値検出手法ごとに「外れている」程度を独自に定量化しています。そのため、手法ごとに検出される外れ値が異なる現象が生じることになります。
外れ値を統一的に定量化できない以上、この出題「100個の玉の中から貴重な石を一つだけ選び出せ！ 」の答えを一つに決めることはできなくなります。
たとえば、naoya_tさんも代表的な外れ値検出法として言及されている１クラスサポートベクターマシンでこの課題を実行すると、１００番目のデータが外れ値となる可能性があります（SVMのパラメータに依存しますが）。k-meansではまた異なった結果ですし、本記事で取り上げていない別の手法ではまた異なる結果が得られるかもしれません。
naoya_tさんによる解答を読む限り、LOFで得られた結果である８７番目のデータ以外は間違いとされるようです。もしそうであれば、たとえば、１クラスサポートベクターマシンの結果から１００番目を外れ値と答えた場合、なぜこの答えが間違いなのかきちんと説明できる必要があります。しかし、そもそもLOFと１クラスサポートベクターマシンでは外れ値の定量化方法が異なっているため、結論を出すことは難しと思われます。
また、どうでもいいことですが、外れ値に相当する玉が本当に貴重な石であるとは限らないことです。単に普通の玉の外れ値である可能性も捨てきれないのですが、この課題では気にする必要はなさそうです。
今回の課題「100個の玉の中から貴重な石を一つだけ選び出せ！ 」を試験問題の一種と捉えると「悪問」（正解が複数あるため）と言わざるを得ないですが、現実的な（場面で遭遇する）課題として見ると「良問」だと思います。現実世界の問題では答えが一つに決まらないことも多々あるという意味で。
というわけでこの外れ値検出課題を解き終わりました。
naoya_tさんは「おわりに」で
「金とキノコ」3部作、いかがでしたでしょうか。やさしすぎると思われた方もいらしたかもしれません。
と書かれていますが、実はやさしくはないのではないかと思いました。まあ、個人的にはとても面白い課題だったのですが。
「金塊か、キノコ料理か」（外れ値検出問題）を解く［LOF(local outlier factor)］
LOF(local outlier factor)とは密度ベースの外れ値検出法です。ある点のまわりの密度がほかの点と比べて小さければ小さいほど、LOFの値は大きくなります。したがって、LOFの最も大きいデータを外れ値すればいいことになります。
LOFアルゴリズムについては後述のアルゴリズム概要に書きましたので興味ある方は参考にしてください。
早速結果を見てみましょう。図から明らかなように、８７番目のLOFが突出して大きいので、これが外れ値となります。
library(DMwR) max.k <- 50 gold.lof <- matrix(0,max.k,100) for(i in 2:max.k){ gold.lof[i,] <- lofactor(gold, k=i) # ここでLOFを計算しています。kはパラメータです。 }
なお、このグラフはLOFのパラメータ$k$が１０～２０のときの結果です（パラメータ$k$については後述のアルゴリズム編を参照してください）。
ところで、naoya_tさんは「「機械学習基礎」簡単な問題を解いて理解しよう！後篇」で、LOFについて次のように説明されています。
LOF（Local Outlier Factor）とは
それぞれのデータ点について、その点からk番目に近い点までの距離を求め、その点のスコアとする
異常値は近傍にデータ点が少ないため、スコアが大きくなる
こうして得られたスコアが基準値を超えるものを検出する手法です。
naoya_tさんは、LOFのことを距離によるスコアをもとに外れ値を求めるアルゴリズムと書かれているようですが、実際は異なります。LOFは距離をベースとした外れ値検出法ではなく、密度をベースとした外れ値検出法です。事実、距離ベースの方法では外れ値をうまく検出できないという問題意識から、LOFは密度ベースの方法を採用しています。
さらに言えば、この文章からはLOFアルゴリズムには外れ値の判断基準値が内在する、または客観的な基準値が存在するようにも読めますが、そのような基準値はなく、外れ値の判定基準はユーザーが決める必要があります。ただし、後述のアルゴリズムを見ればわかるように、あるデータのLOFが１よりかなり大きければ、そのデータが外れ値だろう、とは言えますが、明確な基準値があるわけではありません。
アルゴリズム概要
LOFのアルゴリズムについて簡単に説明します。詳細は参考文献を参照してください。
まず次の記号を導入します。
$D$：データセット
$k$：任意の正数
$d(p,o)$：データ$p \in D$と$o \in D$の距離
なお、この$D$はRの関数lofactor(data, k)の引数dataに、$k$はkに対応します。
定義１
次の関係を満たすデータ$o \in D$と$p$の距離を、$p$の$k$距離（$k-distance(p)$と記述）といいます。
(i)少なくとも$k$個のデータ$o' \in D-\{p\}$に対して、$d(p,o') \leq d(p,o)$が成り立つ
(ii)高々$k-1$個のデータ$o' \in D-\{p\}$に対して、$d(p,o') < d(p,o)$が成り立つ
定義２
$p$の$k$距離が与えられたとき、$p$の$k$距離近傍[$k$-distance neighborhood of $p$]（$N_{k-distance(p)}(p)$または$N_k(p)$と記述）を次で定義します。
$$
N_{k-distance(p)}(p)=\{q \in D-\{p\}|d(p,q) \leq k-distance(p)\}
$$
すなわち、$N_{k-distance(p)}(p)$は$p$の$k$最近傍の含まれるデータのことです。
定義３
データ$p$の$o$に関する到達可能距離[reachability distance of $p$ with respect to $o$]（$reach-dist_{k}(p,o)$と記述）を次の式で定義します。
$$
reach-dist_{k}(p,o)=\max\{k-distance(o), d(p,o)\}
$$
データ$p$が$o$から離れているときには、両データ間の到達可能距離は単純に両データ間の距離となります。一方、両データが十分近いときには、距離の代わりに$o$の$k$距離が到達可能距離となります。
このようにすると、データ$p$が$o$に近くにある場合の距離$d(p,o)$の変動を緩和することができます。
定義４
$p$の局所到達可能密度[local reachability density of $p$]（$lrd_{k}(p)$と記述）を次の式で定義します。
$$
lrd_{k}(p)=\frac{|N_k(p)|}{\sum_{o \in N_k(p)}reach-dist_{k}(p,o)}
$$
すなわち$p$の局所到達可能密度とは、$p$の$k$最近傍あるデータの到達可能距離の平均の逆数です。
以上で準備が整いましたので、LOF(local outlier factor)を定義しましょう。
定義５
$p$の局所外れ係数［local outlier factor of $p$]（$LOF_{k}(p)$と記述）を次の式で定義します。
$$
LOF_{k}(p)=\frac{\sum_{o \in N_k(p)}\dfrac{lrd_{k}(o)}{lrd_{k}(p)}}{|N_k(p)|}
$$
$LOF(p)$は、$p$の局所到達可能密度と$p$の$k$近傍の局所到達可能密度との平均で、$p$の外れ度合いを表しています。
$p$の局所到達可能密度が小さく、そして、$p$の$k$近傍の局所到達可能密度が高くなるにつれて、$LOF(p)$は大きくなります。
LOFの性質として、かたまっている（密集している）データのLOFは１に近いことが知られています。すなわち、LOFが１に近いデータは外れ値の可能性が低く、１からかなり離れているデータでは外れ値の可能性が高くなります。
実際、課題データのLOFのグラフを見ると、多くのデータのLOFがほぼ１であることがわかります。
さて、定義５からわかるように、LOFにはパラメータ$k$（正の整数）があり、その値はユーザーが指定します。
$k$にはどのような値を指定すればよいのでしょうか。
$k$を変化させると、LOF値が変動します。そこでLOF開発者たちはLOF値の標準偏差の変動が安定し始める時点を$k$の最小値としています。
そこで、開発者と同様な方法で今回の課題データに対するパラメータ$k$の範囲を検討してみました。下図は、LOFの標準偏差のグラフとLOFの最小値、最大値、平均値をプロットしたグラフです（標準化した課題データの場合）。
標準偏差のグラフから判断して、$k=$14～20が妥当と思われます。
一方、naoya_tさんは元の課題データを標準化されていませんので、その場合の$k$の範囲がどうなるか検討しましょう。次の標準偏差のグラフ（Indexは２から始まることに注意してください）から、およそ$k=$8～15くらいでしょうか。
naoya_tさんは「ここでは k=3 としていますが、本設問ではk=2～10ぐらいの範囲なら同じ玉がトップに来るようなデータになっています。」と書かれていますので、上述の手順ではなくエイヤッと$k$を決められたのかもしれません。
ところで、論文では$k$の推奨値として、１０～２０を推奨しています。
標準化していないデータに対して、このガイドラインに従った$k$を指定するとどうなるでしょうか。次の図は、$k$（横軸）を変化させたときの外れ値の番号（縦軸）のグラフです。
$k$が１３以上になると外れ値が正解の８７番から１番に変わってしまいます。回答者の中には運悪く、$k=13$以上でLOFを計算してしまい、１番を外れ値とした人もいるかもしれません。
一方、標準化したデータはではどうなるでしょうか。$k$が２８以上になると同じく１番が外れ値となりますが、標準偏差による$k$の範囲（14～20）には含まれまれていませんので解答を間違うことはないですね。
参考文献
Breunig, M. M.; Kriegel, H. -P.; Ng, R. T.; Sander, J. (2000). "LOF: Identifying Density-based Local Outliers". ACM SIGMOD Record 29: 93
「金塊か、キノコ料理か」（外れ値検出問題）を解く［ランダムフォレスト］
一般的に、ランダムフォレストは分類や回帰問題に用いられますが、実はデータ間の近接度も求めることができます。この近接度から外れ度（後述）を計算できるので、この値が大きいデータを外れ値とみなすことができます。
Rスクリプトとその結果は次のとおりです。８７番目のデータが外れ値になります。なお、実行ごとに結果が異なるため次のスクリプトでは２０回実行しています。
下の図は外れ度のグラフ（一例）です。
library(randomForest) outliers <- c() for(i in 1:20){ gold.rf <- randomForest(gold, proximity=TRUE) outliers <- c(outliers,order(outlier(gold.rf),decreasing=TRUE)[1:3]) } table(outliers) outliers 16 17 18 73 79 87 1 1 15 17 6 20
アルゴリズム概要
$prox(n,k)$をデータ$n$と$k$の近接度、$cl(k)$をデータ$k$が属するクラス、$N$をデータ数とします。このときデータ$n$の外れ度（outlier measure）は次の式で定義されます（詳細は参考文献参照）。
$$
\frac{O(n)-\text{median}(O(n))}{||O(n)-\text{median}(O(n))||}
$$
ここで、$O(n)$は
$$
O(n)=\frac{N}{\sum_{cl(k)=j}prox(n, k)^2}
$$
です。
この外れ度が大きいデータが外れ値となります。
参考文献：
