Coursera Machine LearningをPythonで実装 - [Week7]サポートベクターマシン（SVM）

前回に引き続き、今回はサポートベクターマシン（SVM）をPythonで実装していきます。

SVMの水面下でやっているアルゴリズムは、理論的に相当最適化されており複雑です。このコースの演習問題でも既に完成済みのSVMの関数を使ってやったので、今回はScikit-learnの組み込み関数にまかせて、自分で実装することはしません。Andrew Ng先生も「SVMを自分で実装することはおすすめしない」と言っていたのでこれで良いでしょう。途中まで実装してみましたがあまりに長くて投げました。
アルゴリズムを[Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%B5%E3%83%9D%E3%83%BC%E3%83%88%E3%83%99%E3%82%AF%E3%82%BF%E3%83%BC%E3%83%9E%E3%82%B7%E3%83%B3)で覗いてみましたが、クーン・タッカーとか久しぶりに見ました。なにやってるんだかは詳しくは理解できませんでしたが、コースで紹介されている程度の「だいたいこんなことやっているのだろう」ぐらいの理解で十分だと思います。綺麗に理解できたら美しさに気づけると思いますが、そこまで頭が足りませんでした。

**追記** gistからソースファイルをダウンロードできます
2次元のSVM　https://gist.github.com/koshian2/9295b443bfb9cb709cdc6889c839ef1d
スパムメールの分類　https://gist.github.com/koshian2/f4ab592c0807e1f142b720e860027898

## これまでの目次
* [Coursera Machine LearningをPythonで実装 - [Week2]単回帰分析、重回帰分析](https://qiita.com/koshian2/items/894a8ca6f2d5ec9ab70a)
* [Coursera Machine LearningをPythonで実装 - [Week3]ロジスティック回帰](https://qiita.com/koshian2/items/eeef29f28dbde1803e6c)
* [Coursera Machine LearningをPythonで実装 - [Week4]ニューラルネットワーク(1)](https://qiita.com/koshian2/items/032847b6b7cd1315e8b8)
* [Coursera Machine LearningをPythonで実装 - [Week5]ニューラルネットワーク(2)](https://qiita.com/koshian2/items/73854bc1f255a76e7ea3)
* [Coursera Machine LearningをPythonで実装 - [Week6]正則化、Bias vs Variance](https://qiita.com/koshian2/items/3795c0d1bf2d86364f07)

##要点まとめ
SVMのカーネル選択の基準と、正則化パラメーターについてメモ代わりにまとめておきます。

###カーネル選択の基準
ロジスティック回帰とSVMどっちを使うかという話。講義より。

|   変数の数（n）   |   データの数（m）   |                             使うモデル、カーネル                            |
|:-----------------:|:-------------------:|:---------------------------------------------------------------------------:|
| 多い（1万など）   | 少ない（10～1000）  | ロジスティック回帰、カーネルなしSVM                                         |
| 少ない（1～1000） | ほどほど（10～1万） | ガウシアンカーネルのSVM（**一番SVMが有効に機能する領域**）                      |
| 少ない（1～1000） | 多い（5万以上）     | 変数を作るか加えるかしてから、ロジスティック回帰またはカーネルなしSVMを使う |

ニューラルネットワークはどの領域でも有効に機能するが、訓練が遅い。また、**SVMは凸最適化なので、局所解を探索していけば最適解にたどり着くことが保証される**とのこと。SVMのコスト関数はminかmax関数（どちらも凸関数）の正の線形和であることから、コスト関数も凸関数であるため。[ロジスティック回帰と同様](https://qiita.com/koshian2/items/90ac367089b5617519b2)かと思われます。

###正則化のパラメーター
ガウシアンカーネルで使うパラメーターのメモ。後述の内容の再掲です。

|                傾向                |          有効な場面          | 正則化のC | RBFのγ |
|:----------------------------------:|:----------------------------:|:---------:|:------:|
| バイアスが低くなり、分散が高くなる | アンダーフィッティングを解消 | 大        | 大     |
| バイアスが高くなり、分散が低くなる | オーバーフィッティングを解消 | 小        | 小     |

ただし$C=\frac{1}{\lambda}, \gamma=\frac{1}{2\sigma^2}$

**数字を大きくすればオーバーフィットする、数字を小さくすればアンダーフィットする**

## 2次元のSVM
###1.データのロード
今回もデータが多いのでハードコーディングはしていません。同一ディレクトリにOctaveのmatファイルがあることを前提とします。

```python
import numpy as np
import matplotlib.pyplot as plt

from scipy.io import loadmat
from sklearn.svm import SVC

# ex6data1をロード
def load_data(dataname):
    data = loadmat(dataname)
    if dataname == "ex6data3":
        return np.array(data['X']), np.array(data['Xval']), np.ravel(np.array(data['y'])), np.ravel(np.array(data['yval']))
    else:
        return np.array(data['X']), np.ravel(np.array(data['y']))

X, y = load_data("ex6data1")
```

ex6dataの1,2には```X,y```が3には```X,Xval,y,yval```が入っています。3のは交差検証データも含みます。

###2.data1のプロット
```python
def plot_data(X, y):
    pos = y == 1
    neg = y == 0
    plt.plot(X[pos, 0], X[pos, 1], "k+", linewidth=1, markersize=7)
    plt.plot(X[neg, 0], X[neg, 1], "ko", color="y", markersize=7)

# X,yをプロット
plot_data(X, y)
plt.show()
```
![ml_wk7_1.png](https://qiita-image-store.s3.amazonaws.com/0/232088/a500d33d-0761-09c0-d795-240eaa6d6c2f.png)
だいたいこんなデータですよと。

###3.data1を線形カーネルのSVMで分類
正則化のパラメーターC（＝ロジスティック回帰のλの逆数）を変化させながら、決定境界をプロットしてみます。カーネル関数は線形関数です。「カーネルを入れないよ」と言ったときは線形カーネルを使うのと同じです。

```python
## 線形カーネル
# C=の値を変えて決定境界をプロット
fig = plt.figure(figsize = (8, 8))
fig.subplots_adjust(hspace = 0.2, wspace = 0.2)
for i, c in enumerate([1, 10, 100, 1000]):
    clf = SVC(C = c, tol=1e-3, kernel="linear")
    clf.fit(X, y)
    w = clf.coef_[0]
    b = clf.intercept_[0]
    xp = np.linspace(np.min(X[:, 0]), np.max(X[:, 1]), 100)
    yp = -(w[0] * xp + b) / w[1]
    # 決定境界をプロット
    ax = fig.add_subplot(2, 2, i+1)
    ax.plot(X[y == 1, 0], X[y == 1, 1], "k+", linewidth=1, markersize=7)
    ax.plot(X[y == 0, 0], X[y == 0, 1], "ko", color="y", markersize=7)
    ax.plot(xp, yp, color="b")
    ax.set_title(f"Linear SVM C = {c}")
plt.show()
```

SVC（Support Vector Classifier）関数はいろいろなカーネルが使えるので、**kernelパラメーターで指定するのを忘れないようにしましょう。デフォルトではガウシアンカーネル（rbf）になります**。

![ml_wk7_2.png](https://qiita-image-store.s3.amazonaws.com/0/232088/7c1bf28c-bb3a-6275-1122-06b06789331a.png)

Cの値を変えていくとだんだん決定境界が傾いて、”寄っていく”のが確認できます。Cはλの逆数なので、Cを小さくするとアンダーフィット、Cを大きくするとオーバーフィットするようになります。ここは交差検証から決めればよいです（後ほどやります）。

###4.ガウシアンカーネル（RBF）
ただのガウス関数です。講義では次のように示されていました。

$$K_{gaussian}(x^{(i)}, x^{(j)}) = \exp\bigl(-\frac{||x^{(i)}-x^{(j)}||^2}{2\sigma^2}\bigr) $$

また$\gamma=\frac{1}{2\sigma^2}$とおいて、

$$K_{gaussian}(x^{(i)}, x^{(j)}) = \exp\bigl(-\gamma||x^{(i)}-x^{(j)}||^2) $$

と書くこともあります。Scikit-learnのSVCでkernel="rbf"パラメーターのγもこの定義です。いくら背景の実装を知らなくてもいいとは言っても、変数変換をしてあげないといけないのでここは知っておく必要があります。

ブラックボックスのままなのはよくないので、カーネル関数を自分で実装して想定通りの値が出たか確認してみます。実際はこのカーネル関数はSVC側で用意してくれるので、特に実装する必要はありません。

```python
##ガウシアンカーネル
# ガウシアンカーネルのテスト（組み込みであるので特に導入する必要はない）
def gaussian_kernel(x1, s2, sigma):
    sim = np.exp(-np.linalg.norm(np.ravel(x1) - np.ravel(x2)) ** 2 / (2 * sigma ** 2))
    return sim
# sigma=2ではテスト用の値は0.324652になる
x1 = np.array([1, 2, 1])
x2 = np.array([0, 4, -1])
sigma = 2
# テスト用の値で類似度計算
sim = gaussian_kernel(x1, x2, sigma)
print("Gaussian Kernel between x1 = [1; 2; 1], x2 = [0; 4; -1], sigma =", sigma)
print("\t", sim)
print("(for sigma = 2, this value should be about 0.324652)")
print()
```

```:出力
Gaussian Kernel between x1 = [1; 2; 1], x2 = [0; 4; -1], sigma = 2
         0.324652467358
(for sigma = 2, this value should be about 0.324652)
```

σの値が大きくなると、カーネル関数がよりゆったりとした幅になりバイアスが大きくなります。つまり、オーバーフィッティングを解消する際に有効です。σが小さくなると、幅が狭くなりモデルの分散が大きくなり、アンダーフィッティングを解消します。γを使う場合は分数になるので、順序関係が逆転します。

頭がこんがらがってくるので、正則化のパラメーターCも含めて整理します。

|                傾向                |          有効な場面          | 正則化のC | RBFのγ | （←のσの場合） |
|:----------------------------------:|:----------------------------:|:---------:|:------:|:--------------:|
| バイアスが低くなり、分散が高くなる | アンダーフィッティングを解消 | 大        | 大     | 小             |
| バイアスが高くなり、分散が低くなる | オーバーフィッティングを解消 | 小        | 小     | 大             |

ここで、$C=\frac{1}{\lambda}$、λは回帰分析で出てきた正則化のパラメーター、$\gamma=\frac{1}{2\sigma^2}$です。SVCのRBFではσではなくγを使うので、Cも含めて**数字を大きくすればオーバーフィットする、数字を小さくすればアンダーフィットする**と覚えておけば大丈夫そうです。

###5.data2をロード
```python
##データ2をロード
X, y = load_data("ex6data2")

# X, yをプロット
plot_data(X, y)
plt.show()
```

このようなデータで線形カーネルでは綺麗に分類するのが難しいケースです。ガウシアンカーネルを使うと非常に綺麗に分類できます。

![ml_wk7_3.png](https://qiita-image-store.s3.amazonaws.com/0/232088/93964c5e-30bb-bdae-71dc-12e08515c1d4.png)

###6.data2をガウシアン（RBF）カーネルのSVMで分類
ガウシアンカーネルではSVCのkernel引数がrbfになっただけで（デフォルトがrbfなので省略してもOK）特に難しい所はありません。**predict_probaで確率を予測したい場合は、SVCの引数のprobabilityをTrueに**しましょう。これを入れると遅くなりますが、デフォルトのFalseのままpredict_probaを叩くと怒られます（要は高速化のためにあえて切っているっぽい）。

```python
## RBFカーネルでSVMを訓練
C, sigma = 1, 0.1
# γ=1/(2σ^2)
clf = SVC(C = C, gamma=1/(2*sigma**2), kernel="rbf", tol=1e-3, probability=True)
clf.fit(X, y)
# 決定境界のプロット
def visualize_boundary(X, y, model):
    # 訓練データのプロット
    plot_data(X, y)
    # グリッドを作り予測値を計算
    x1plot = np.linspace(np.min(X[:, 0]), np.max(X[:, 0]), 100)
    x2plot = np.linspace(np.min(X[:, 1]), np.max(X[:, 1]), 100)
    X1, X2 = np.meshgrid(x1plot, x2plot)
    vals = np.zeros(X1.shape)
    for i in range(X1.shape[1]):
        this_X = np.c_[X1[:, i], X2[:, i]]
        vals[:, i] = model.predict_proba(this_X)[:, 1]
    # 等高線プロット
    plt.contour(X1, X2, vals, levels=[0.5])

visualize_boundary(X, y, clf)
plt.show()
```

predict_probaでは列方向にy==0である確率、y==1である確率を計算しているので、y==1である確率のみを選択して0.5をしきい値として決定境界をプロットします。

綺麗で美しい決定境界が描けました。ロジスティック回帰で多項式変数を大量にぶちこんでオーバーフィッティングしたときのギザギザはありません。なるほど、これは頼もしいですね。

![ml_wk7_4.png](https://qiita-image-store.s3.amazonaws.com/0/232088/02ede975-8d3c-9066-cd9d-e6781cfdbf4e.png)

###7.data3のロード
data2ではCとσをアドホックに与えましたが、そもそもこの2つはハイパーパラメータなので、交差検証データからグリッドサーチして誤差が最小化されるような値を与えるのが適切です。data2よりデータを少なくしたものがdata3です。

```python
## データ3をロード
#～valは交差検証データ
X, Xval, y, yval = load_data("ex6data3")
plot_data(X, y)
plt.show()
```

X,yが訓練データ、Xval,yvalは交差検証データとなります。

![ml_wk7_5.png](https://qiita-image-store.s3.amazonaws.com/0/232088/00dbbc7f-a618-828f-3b6b-5851ee59a812.png)

###8.data3をRBFのSVMで分類＋ハイパーパラメータのグリッドサーチ
グリッドサーチはわざわざこのようなfor文を書かなくても、[Scikit-learnの組み込みでできるようです](https://qiita.com/SE96UoC5AfUt7uY/items/c81f7cea72a44a7bfd3a)。ただの総当たり法なので、とりあえずfor文で書いてみました。

```python
## RBFカーネルでSVMを訓練
# ハイパーパラメータの調整
def dataset3_params(X, y, Xval, yval):
    C, gamma = 1, 1/(2*0.3**2)
    error = np.inf
    params = np.array([0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30])

    for ci in params:
        for si in params:
            model = SVC(C=ci, gamma=1/(2*si**2))
            model.fit(X, y)
            pred = model.predict(Xval)
            error_tmp = np.mean(pred != yval)
            if error_tmp < error:
                error, C, sigma = error_tmp, ci, si
    return C, sigma
# ハイパーパラメータを調整する
C, sigma = dataset3_params(X, y, Xval, yval)
print(f"C = {C}, sigma = {sigma} is selected from cross validation set.")
# 選択したハイパーパラメータをもとに訓練
# visualize_boundaryの中でpredict_probaを使っているので、probability=Trueとする必要がある
model = SVC(C=C, gamma=1/(2*sigma**2), probability=True)
model.fit(X, y)
visualize_boundary(X, y, model)
plt.show()
```

やっていることは単純で、Cとσ（γ）を総当たりで計算して、誤差が最小となる組み合わせを返しているだけです。errorの値を最初無限大とおいて、小さくなれば更新…と実装すればすぐできますね。こんなコードを書かなくても[組み込み](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html)を使えばもっと明瞭にできますし、異なるカーネル間のサーチもできるようです。

```:出力
C = 1.0, sigma = 0.1 is selected from cross validation set.
```
C=1.0, σ=0.1が選ばれました。このときの決定境界は次の通りです。

![ml_wk7_6.png](https://qiita-image-store.s3.amazonaws.com/0/232088/912be21f-aab4-ad8f-3667-140a58fc1163.png)

いい感じにできてそう。

## スパムメールの分類
演習の2つ目でより実践的なものとして、スパムメールをSVMによって分類する問題がありました。前処理までやっているので「こんな感じでやるのかー」と参考になります。

###1.おおまかな流れ
SVMの問題に落とし込むために、メールの本文を次のようなベクトルに変換します。

1. 事前処理で固有名詞やHTMLタグ、電話番号など、あまり必要なさそうな情報を落とす
2. 単語単位で分割（日本語みたいに形態素解析しなくていいのが楽すぎる）
3. ステミング（[IBMのページより参考](https://www.ibm.com/support/knowledgecenter/ja/SS8NLW_11.0.1/com.ibm.swg.im.infosphere.dataexpl.engine.doc/c_correcting_stemming_errors.html)）という、語幹を抽出する操作を行う。日本語で言う所の活用形を統一するみたいな感じだと思われる。
4. 単語リスト（これは事前に用意されたもので、よく出てくる単語を集めたもの）から、本文に含まれる単語のインデックスを求める。例えば、ourが500番で、priceが700番だったらこの500と700の数字をメモする。
5. 4で求めた数字をもとに、単語リストの長さ分の0ベクトルを用意し、500と700の要素を1にする。これを多数用意したものが訓練データとなる。

ほとんど前処理が重要な気がする。ここらへんは分析する人がやりやすい言語でやってもよさそう。

###2.単語リストの作成
vocab.txtというテキストファイルが、

```
10	accept
11	access
12	accord
13	account
```

こんな感じでタブ区切り＋改行で並んでいたとします。まずは単語リストを作ります。ちょっとオブジェクト指向風に。

```python
import re
import numpy as np

from nltk.stem import PorterStemmer
from scipy.io import loadmat
from sklearn.svm import SVC

## 単語リストを作る
# 単語用のクラス定義
class Vocab:
    def __init__(self, str_line):
        spl = str_line.split()
        self.index = int(spl[0])
        self.word = spl[1]
# 単語リスト
vocab_list = []
# ファイル読み込み
with open("vocab.txt", "r") as fp:
    vocab_data = fp.readlines()
# 単語のインスタンスを作成
for v in vocab_data:
    vocab_list.append(Vocab(v))
```

###3.メールの前処理
```python
## メールデータの前処理
## →メールデータを解析しやすい形に変換して、単語リストに入っているインデックスを返す
def process_email(email_contents):
    # 単語リストのグローバル変数を使う
    global vocab_list
    # 返り値の初期化
    word_indices = []

    ## 前処理
    # 全て小文字に変換
    contents = email_contents.lower()
    # HTMLタグをスペースに変換
    contents = re.sub(r"<[^<>]+>", " ", contents)
    # （電話）番号をnumberという単語に変換
    contents = re.sub(r"[0-9]+", "number", contents)
    # URLをhttpaddrという単語に変換
    contents = re.sub(r"(http|https)://[^\s]*", "httpaddr", contents)
    # メールアドレスをemailaddrという単語に変換（真ん中の@を探す）
    contents = re.sub(r"[^\s]+@[^\s]+", "emailaddr", contents)
    # $マークをdollarという単語で統一する
    contents = re.sub(r"[$]+", "dollar", contents)

    # 単語単位に分割する
    contents_words = re.split("[\'\s@$\/#.-:&*+=\[\]?!(){},\">_<;%]", contents)
    # stemmerのインスタンスを作る
    stemmer = PorterStemmer()
    # ヘッダー
    print("\n==== Processed Email ====\n")
    # 文字数
    l = 0

    ## 単語単位でサーチ
    for str in contents_words:
        # アルファベット以外を取り除く
        str = re.sub(r"[^a-zA-Z0-9]", "", str)
        # PorterStemmerで語幹を取り出す（stemming）
        # nltk.stemに用意されており、Anacondaでインストールすると入っているはず
        str = stemmer.stem(str.strip())

        # 空白文字か短すぎる文字はスルー
        if len(str) < 1: continue

        # ルックアップする
        query = list(filter(lambda vocab: vocab.word == str, vocab_list))
        if len(query) == 1:
            word_indices.append(query[0].index)
        # 適度に改行
        if (l + len(str) + 1) > 78:
            print()
            l = 0
        # 画面に表示
        print(str, end=" ")
        l += len(str) + 1
    # フッター
    print("\n\n=========================\n")
    # 返り値
    return word_indices
```

グローバル変数を使うよ宣言をして、まずは正規表現でがーっと置き換えています。正規表現はOctaveもPythonも同じなので、移植が簡単。

また正規表現を使って単語単位に分割し、```PorterStemmer()```というのがステミングのインスタンスです。```nltk```というパッケージに入ってたもので、Anaconda経由でインストールしていたら多分インストールされているはずです。すぐこういうのが出てくるのがPython強い。

あとは特に言うまでもなくルックアップしてるだけです。

###4.前処理の確認
サンプル（emailSample1.txt）を読んで、前処理が上手く行っているか確認します。ちなみにこんな内容。

```:emailSample1.txt
> Anyone knows how much it costs to host a web portal ?
>
Well, it depends on how many visitors you're expecting.
This can be anywhere from less than 10 bucks a month to a couple of $100. 
You should checkout http://www.rackspace.com/ or perhaps Amazon EC2 
if youre running something big..

To unsubscribe yourself from this mailing list, send an email to:
groupname-unsubscribe@egroups.com
```

```python
## メールの前処理
print("Preprocessing sample email (emailSample1.txt)")
# 特徴を抽出
with open("emailSample1.txt") as fp:
    file_contents = fp.read()
word_indices = process_email(file_contents)
# 画面に表示
print("Word Indices: ")
print(word_indices)
print()
```

```:出力
Preprocessing sample email (emailSample1.txt)

==== Processed Email ====

anyon know how much it cost to host a web portal well it depend on how mani
visitor you re expect thi can be anywher from less than number buck a month
to a coupl of dollarnumb you should checkout httpaddr or perhap amazon ecnumb
if your run someth big to unsubscrib yourself from thi mail list send an
email to emailaddr

=========================

Word Indices:
[86, 916, 794, 1077, 883, 370, 1699, 790, 1822, 1831, 883, 431, 1171, 794, 1002,
 1893, 1364, 592, 1676, 238, 162, 89, 688, 945, 1663, 1120, 1062, 1699, 375, 116
2, 479, 1893, 1510, 799, 1182, 1237, 810, 1895, 1440, 1547, 181, 1699, 1758, 189
6, 688, 1676, 992, 961, 1477, 71, 530, 1699, 531]
```

原文→前処理→行列変換としていくうちに、どんどんコンピューターにとって解釈しやすい形になっていくのがわかります。モデルに落とし込むというのはこういうことらしい。

###5.インデックス→ブール値のベクトルにマッピング
インデックスのベクトルのままではデータごとに要素数が異なり、うまく行列計算できないので、単語が出現した、出現しないの0or1の要素数一定のベクトルにマッピングします。ここでの要素数は単語リストの数と一致します。

```python
# インデックスを変数にマッピング
def email_feature(word_indices):
    n = 1899
    x = np.zeros((n, 1))
    for i in word_indices:
        x[i-1] = 1
    return x

# 変数のベクトル
features = email_feature(word_indices)
print("Length of feature vector:", len(features))
print("Number of non-zero entries:", np.sum(features > 0))
print()
```

1899次元、横がおおよそ1900次元のベクトルのSVM最適化問題に落とし込むことができました。あとは計算するだけです。

```
Length of feature vector: 1899
Number of non-zero entries: 45
```

###6.SVMを訓練データで計算
既に大量の（訓練データが4000件、テストデータが1000件）前処理済みの行列がspamTrain.mat、spamTest.matに用意されています。まずは訓練データを読み込んでフィットさせます。

```python
## 線形SVMを使ってスパムを分類する
# spamTrain.matには既にベクトル化されたデータが入っている
data = loadmat("spamTrain.mat")
X, y = np.array(data['X']), np.ravel(np.array(data['y']))
C = 0.1
# SVM
print("Training Linear SVM (Spam Classification)")
model = SVC(C=C, kernel="linear")
model.fit(X, y)
# 訓練データの精度
p = model.predict(X)
print("Training Accuracy:", np.mean(p == y) * 100)
```

訓練データではほぼ100%でフィットしていることになります。これで線形カーネルで、計算が一瞬（10秒弱）だから、ガウシアンカーネルを使ったりすればまだまだ余力がありそうなのが素晴らしい。それだけSVMのアルゴリズムが洗練されたものなのでしょう。

```:出力
Training Linear SVM (Spam Classification)
Training Accuracy: 99.825
```

ついでに寄与率の高いトップ15の単語をリストアップしてみます。これらのワードが多いほどスパムメールとして判定されやすいようです。

```python
# 係数の高いトップ15の単語
topword_indices = np.argsort(model.coef_[0])[::-1][:15]
for i in topword_indices:
    print(vocab_list[i].word, model.coef_[0,i])
print()
```

ourはいろんな英文で使うからまだしも、click, guarantee, click, visit, dollar, priceあたりは怪しいですね。リンク誘導、性能や返金保証、料金の話、ここらへんがスパムと判定されやすいのかもしれない。

```:出力
our 0.500613736175
click 0.465916390689
remov 0.422869117061
guarante 0.383621601794
visit 0.367710398246
basenumb 0.345064097946
dollar 0.323632035796
will 0.269724106037
price 0.267297714618
pleas 0.2611688867
most 0.257298197952
nbsp 0.25394145516
lo 0.253466524314
ga 0.248296990456
hour 0.246404357832
```


###7.テストデータの精度
本来ほしいのでは訓練データの精度ではなく、テストデータの精度なのでそれも計算します。

```python
## テストデータの精度
# spamTest.matには既にデータがある
data = loadmat("spamTest.mat")
Xtest, ytest = np.array(data['Xtest']), np.ravel(np.array(data['ytest']))
# テストデータで予測
p = model.predict(Xtest)
print("Test Accuracy:", np.mean(p == ytest) * 100)
```

テストデータの99%近く当てられるらしい。すごい。

```:出力
Test Accuracy: 98.9
```

###8.任意のスパムメールの判定
スパム判定器ができたので、任意の（といっても既に用意されているスパムメールのサンプルですが）メールを判定できることを確認します。判定するのはこれ。

```:spamSample2.txt
Best Buy Viagra Generic Online

Viagra 100mg x 60 Pills $125, Free Pills & Reorder Discount, Top Selling 100% Quality & Satisfaction guaranteed!

We accept VISA, Master & E-Check Payments, 90000+ Satisfied Customers!
http://medphysitcstech.ru
```

バイアグラの宣伝です。どう見てもスパムです本当にありがとうございました、なやつ。

```python
## 任意のデータ
with open("spamSample2.txt") as fp:
    file_contents = fp.read()
word_indices = process_email(file_contents)
features = email_feature(word_indices)
p = model.predict(np.array(features).reshape(1, -1))
print("Spam Classification:", p)
print("(1 indicates spam, 0 indicates not spam)")
```

```:出力
Spam Classification: [1]
(1 indicates spam, 0 indicates not spam)
```

1として判定されたからスパムだよ、とのこと。そりゃのあの文面見ればそうですよね。じゃあちょっと微妙そうなやつをやってみます。ありました。

```:myEmail.txt

Quick Answers about our New Droplet Plans & Pricing

Quick answers about Bandwidth Usage and Billing

Thank you for submitting your ticket. Here’s what we received:

Hello, 
My account was locked suddenly. 
Can you tell me why it was locked? 

>Account Locked 
>Your account has been flagged for abuse. If you feel that our system has made an error please open a support ticket.

Thanks so much,

DigitalOcean Support

ref:*****:ref
Copyright c 2017 DigitalOcean
Floor 10, 101 Avenue of the Americas, New York, NY, 10013
All rights reserved.
```

リファラを一部伏せています。これは、DigitalOceanというサーバー会社に問い合わせたときに（内容は、突然アカウントがロックされたからなぜロックされたのか教えて欲しいというものです）、コピーとして発行されたメールで、**Yahooメールでは迷惑メールとして分類されていました**。本当に迷惑メールとして判定される要素があるのか確認してみます。

```:spamSample2.txtのスパムテスト
Spam test on spamSample2.txt ...

==== Processed Email ====

best buy viagra gener onlin viagra numbermg x number pill dollarnumb free
pill reorder discount top sell number qualiti satisfact guarante we accept
visa master echeck payment number satisfi custom httpaddr

=========================

Spam Classification: [1]
(1 indicates spam, 0 indicates not spam)
Estimated probability of non-spam / spam :
[[ 0.06644329  0.93355671]]
```

```myEmail.txtのスパムテスト
Spam test on myEmail.txt ...

==== Processed Email ====

quick answer about our new droplet plan price quick answer about bandwidth
usag and bill thank you for submit your ticket here what we receiv hello my
account wa lock suddenli can you tell me whi it wa lock account lock your
account ha been flag for abus if you feel that our system ha made an error
pleas open a support ticket thank so much digitalocean support ref
numberdonumber**** number****** ref copyright c number digitalocean floor
number number avenu of the america new york ny number all right reserv

=========================

Spam Classification: [1]
(1 indicates spam, 0 indicates not spam)
Estimated probability of non-spam / spam :
[[ 0.44956828  0.55043172]]
```

今回はSVCのProbability=Trueにして確率も計算できるようにしました。**確かにYahooメールの判定はある意味正しくて、このスパム判定機でもスパムであると判定されてしまいました**。ただ確率を見ると、前者はバイアグラの宣伝という誰がどうみてもスパムなやつなので、スパムである確率は高く93.4%であるのに対して、後者はスパム確率55.0%と割と微妙な判定でした。ただ、ここの会社[サポートが割と不親切](https://www.lowendtalk.com/discussion/106692/digitalocean-abuse-account-problem)なので（チケットへの返信は遅い割に[Twitterでぼやくと凄い速度でクソリプ送ってくる](https://twitter.com/koshian2/status/991674341422071810/photo/1)）、Yahooの判定が正しかったのが証明されたのはちょっとニヤリとしました。

以上です。今回スパムメールという初めて実践的な例をやりましたが、前処理がかなり大変そうな印象を受けました。