lazy diary

統計とその周辺

尤度比検定、Wald検定、Score検定の着目する点

この3つ検定統計量が違うだけで結局何が違うんだと思ったのでメモ。

概要

・前提知識

以下、パラメータ\thetaの次元は1次元としておきます。簡単に定式化しておきます。
独立な標本X_1, X_2, \dots, X_n \sim f(x;\theta)としておきます。この標本から得られる最尤推定量(以下、MLEとします)を\hat{\theta}_nとしましょう。
このとき確率密度関数f(x;\theta)がいくつかのある仮定を満たすとき、MLEは漸近正規性を持ちます。すなわち、


\sqrt{n I_1(\theta_{0})} (\hat{\theta}_{n}-\theta_{0}) \to N(0,1)
のように分布収束します。I_n(\theta)はフィッシャー情報量です。nは標本の数を表します。

・尤度比検定

対数尤度に着目した検定です。
帰無仮説と対立仮説 H_0 :\theta \in \Theta_0,  \quad H_1 : \theta \in \Theta_0^c
対数尤度 L(\theta | X)
このとき次の写像を尤度比統計量と言います。


\lambda(x)=\frac{\sup_{\theta \in \Theta_0} L(\theta | X)}{\sup L(\theta | X)}
このとき尤度比統計量の分母はただのsupなので基本的には\theta = \hat{\theta}_nとMLEを代入した値となります。このとき

\lambda(x) \to \chi_1^2
のように分布収束します。この分布を使う検定です。

・Wald検定

パラメータに着目した検定です。
帰無仮説と対立仮説 H_0 :\theta = \theta_0,  \quad H_1 : \theta \neq  \theta_0
このときMLEの漸近正規性とスラツキーの定理から


\sqrt{nI_1(\hat{\theta}_n)}(\hat{\theta}_n - \theta) \to N(0, 1)
のように分布収束します。この分布を使う検定です。

・Score検定

スコア関数(つまり対数尤度の導関数)に着目した検定です。
スコア関数S_n(\theta, X) = dL/d\theta
ある仮定のもとスコア関数の期待値は0で、さらに分散はフィッシャー情報量と等しくなることから


\frac{S_n(\theta, X)-E[S_n(\theta, X)]}{\sqrt{V[S_n(\theta, X)]}} = \frac{S_n(\theta, X)}{\sqrt{nI_1(\theta)}} \to N(0, 1)
右辺の分布収束は中心極限定理によるものです。スコア関数が各標本に関するカタマリの和であるため定理が利用できます。




Score検定だけMLEを必ずしも必要としません(帰無仮説H0で定めた値を用いるだけ)。加えてScore検定はスコア関数の性質を使いMLEの漸近正規性までは使わないので少し仮定が緩いはずです。となるとScore検定を使うのが一番縛りが少ないのかなと思いました。参考文献2によると漸近的には3つの検定どれも等しいらしいが果たして。。。

参考文献

1. 久保川達也 『現代数理統計学の基礎』共立出版
6, 7章が上記に関係する話です。

2. https://socialsciences.mcmaster.ca/jfox/Courses/SPIDA/mle-mini-lecture-notes.pdf
3つの検定の観点の違いがグラフで一目で分かるようになっている。

統計検定1級に合格しました

合格したので感想を。合格した勉強法の1例として参考になれば。
2018年11月下旬に統計検定1級を受験して統計数理と統計応用(理工学を選択)の両方に合格しました。以下、その受験記です。
http://www.toukei-kentei.jp/wp-content/uploads/ToukeiAkari-300x123.png*1



背景

自分のことを少しだけ書きます。数学を専攻している大学生です。とは言っても大学で統計学を学んだことはなく、2018年3月ごろから趣味で勉強を始めました。2018年5月には統計検定2級に合格しています。良ければその受験記も見て下さい。
統計検定2級に合格しました - lazy diary
さて、「2級に合格したなら次は1級でしょ」ということで受験しました。準1級を飛ばしたのは単純に時期が合わなかったからです。


統計検定1級の概要

年1回11月下旬開催の試験です。試験範囲をチェックしましょう。
www.toukei-kentei.jp
試験内容は統計数理と統計応用の2つに分かれています。
この両方に合格すると統計検定1級合格となります。

・統計数理

試験日の午前に行なわれる90分の筆記試験です。60/100点で合格です。
5つの問題から3問を選んで解きます。確率分布に関わる様々な計算がメインです。
確率密度関数の計算、確率母関数を用いた計算、検定論が頻出です。

・統計応用

試験日の午後に行なわれる90分の筆記試験です。60/100点で合格です。
統計応用は「人文科学」「社会科学」「理工学」「医薬生物学」の中から1つの分野を試験申し込み時に選んで解くことになります。
事前に選んだ分野の5つの問題から3問を選んで解きます。



勉強したこと

さて、まず参考書選びです。1級に合格したというブログを調べるとちらほら出てくるのでそこにオススメしてある参考書と、あとは公式サイトにある合格者の声にリストアップされている参考書を書店でチェックしました。
www.toukei-kentei.jp
その中で私が選んだのがこちらの本。

・現代数理統計学の基礎 (共立講座 数学の魅力)

現代数理統計学の基礎 (共立講座 数学の魅力)

現代数理統計学の基礎 (共立講座 数学の魅力)

この本が無ければ合格は無かったです。それくらい良い本でした。 *2
統計検定1級統計数理の範囲となっているが本には掲載されていない事項の一部を著者がWebページに補足資料として載せています。それくらい統計検定に最適化された本です。この本で統計数理のほとんどの試験範囲を抑えることができます。以下、『現代数理統計学の基礎』の目次です。

第1章 確率
第2章 確率分布と期待値
第3章 代表的な確率分布
第4章 多次元確率変数の分布
第5章 標本分布とその近似
第6章 統計的推定
第7章 統計的仮説検定
第8章 統計的区間推定
第9章 線形回帰モデル
第10章 リスク最適性の理論
第11章 計算統計学の方法
第12章 発展的トピック:確率過程
引用:https://www.kyoritsu-pub.co.jp/bookdetail/9784320111660

第1章から第9章までは十分に理解しておくべきでしょう。試験問題の半分くらいは適切に確率密度関数や期待値、分散などの計算ができるかどうかです。この計算がスムーズにできるかどうかが鍵だと思います。
第10章「リスク最適性の理論」は試験範囲の観点では飛ばしても問題ないと思います。
第11章、第12章についても理解しておくべきですが、試験に出ないであろう内容もあるので試験範囲と照らし合わせて必要に応じてとい感じで。

この本をノートにまとめながら読んで、さっそく過去問に取り組みます。2017年の統計数理の問題を解きました。自己採点をしたら合格ライン上の点数だったと思います。それより時間が足りません。明らかに演習不足なのを実感しましたが、同時に解説を読めば理解できるということも分かったので統計数理の勉強はひとまず終了にしました。



次に統計応用の勉強をします。私は「理工学」分野を選択しました。理由は統計数理とさほど範囲が変わらないからです。過去問を見てもらえば分かると思います。「理工学」以外の分野に明るい場合はそれを選べば良いと思います。特にそういうのがないですという場合は、「理工学」分野を選ぶのが無難だと私は思います。
参考書は次の2冊です。

・多変量解析入門――線形から非線形

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

1 はじめに
2 線形回帰モデル
3 非線形回帰モデル
4 ロジスティック回帰モデル
5 モデル評価基準
6 判別分析
7 ベイズ判別
8 サポートベクターマシーン
9 主成分分析
10 クラスター分析
引用:https://www.iwanami.co.jp/book/b265441.html

この本はかなり丁寧に記述されています。統計検定云々は関係なく、多変量解析を学びたい人にはオススメしたい1冊でした。 *3

・日本統計学会公式認定 統計検定1級対応 統計学

日本統計学会公式認定 統計検定1級対応 統計学

日本統計学会公式認定 統計検定1級対応 統計学

  • 作者: 二宮嘉行,大西俊郎,小林景,椎名洋,笛田薫,田中研太郎,岡田謙介,大屋幸輔,廣瀬英雄,折笠秀樹,日本統計学会,竹村彰通,岩崎学
  • 出版社/メーカー: 東京図書
  • 発売日: 2013/04/08
  • メディア: 単行本
  • この商品を含むブログ (5件) を見る

第I部 統計数理
第1章 確率と確率変数
第2章 種々の確率分布
第3章 統計的推定
第4章 仮説検定
第5章 主なデータ解析手法

第II部統計応用
第6章 統計応用共通手法
第7章 人文科学分野キーワード
第8章 社会科学分野キーワード
第9章 理工学分野キーワード
第10章 医学生物学分野キーワード
引用:http://www.tokyo-tosho.co.jp/books/ISBN978-4-489-02150-3.html

次は公式の参考書です。上2冊で学習できていない分散分析や直交表などを補いつつ、これまで学んだ事を一気に復習します。*4



さて、これで統計応用の範囲も全てチェックしました。
私は時間にゆとりがあったのでこの間にもいくつか別の本をやっていますが、これだけで十分だと思います。後はひたすら演習です。時間をきっちり計って過去問を解きまくりましょう。
試験日まで時間がある人はこのあたりで気づくと思うのですが、年1回しか試験がないので過去問が少ないです。もっと演習をしたいという人にはこれが良いかもしれません。
Exam past papers
王立統計学会(RSS)が実施している試験の過去問です。これも統計検定公式サイトに載っている合格者の声で知りました。このうち「Graduate Diploma」の過去問を解いてみるのが良いと思います。難易度は1級よりやや簡単ですが、問題に慣れるのには十分でしょう。実際に似たような問題が1級試験でも出題されています。



受験当日

受験票に写真を貼ることと、電卓、時計を忘れないようにしましょう。
会場に行きます。最初にアンケートがあるので回答しましょう。
さっそく統計数理の試験開始です。どこかで見たことがある問題を見つけ1問サクッと解答。残り2問は最後がダメでしたがそれ以外は一通り解けたと思います。そしてお昼休憩を挟んで、午後から統計応用の試験です。
計算量が多そうな条件付き正規分布の計算問題はスルーして、とにかく分かる問題を落とさずしっかり得点するという気持ちで解いていきました。

12月半ばにWeb合格発表が行われ無事合格しました。



最後に

統計検定という名前がついていますが試されていることの半分は「数式を素早く記述できますか」という点だと思いました。

統計検定1級を通して、いわゆる古典的な統計は一通り学べたと思っています。今はDeep Learningが流行っていますが、確率関数の操作だったりはそれらを理解するための基本なので、その理解度をチェックする上で統計検定を受験してみてはいかがでしょうか。


・更新履歴
2018/12/17. 記事作成
2019/07/22. 一部表現を変更

*1:統計検定の公式キャラクターらしいです http://www.toukei-kentei.jp/info/toukeiakari/

*2: 私は数学畑の人なので苦ではなかったのですが、これは数式が沢山並んでいる本です。と言っても数理統計なのでそうなるのも仕方ないのですが、普段あまりこの手の本を手に取らない方は読むのに少々苦労するかもしれません。

*3: もし「試験日まで時間がない」とか「とにかく近道して受かりたい」といった場合はこの本は飛ばしても問題ないでしょう。これら多変量解析の内容は準1級および1級試験の試験範囲に含まれています。従って試験に出題されてもおかしくないです。実際に過去問を見てみると準1級試験ではこれらの出題事例があります。例えば、2018年準1級試験ではマーク問題(not記述式)として線形判別やSVMが出題されています。一方で1級試験では私の知る限り出題事例が無く、せいぜい線形重回帰くらいだったと思います。

*4: この本は某レビューで酷評されているようですが、私は使い方の問題だと思っていて、この本でイチから学ぼうとすると確かに難しい気がしますが、復習に使うならコンパクトに試験範囲がまとまっているので良い本だと思います。

非線形SVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

ハードマージンSVM、ソフトマージンSVMと続いて今回でSVMは最後です。
今までは全て線形SVMでしたが、今回は非線形SVMとなります。非線形SVMでは特徴空間次元の内積計算が多数必要になり、それを避けるためにカーネル法という手法が用いられます。カーネル法にも種類は様々あるのですが今回は最も一般的なガウスカーネルを用います。データはソフトマージンSVMのときと同じデータです。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class NonLinearSoftMarginSVM():
  
  def __init__(self, X, y, lam, sigma):
    self.X = X
    self.y = y
    self.lam = lam
    self.sigma = sigma
    
  #ガウスカーネル
  def kernel(self, x, y):
    return np.exp(np.dot(x - y, x - y) / (-2 * self.sigma))
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * self.kernel(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      for j in range(num_of_data):
        b = b + alpha[j] * label.iloc[j] * self.kernel(self.X.iloc[j], self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index, alpha]

前回のソフトマージンSVMと比べて変化しているのは、内積計算が全てカーネルによる計算に置き換えられているという点だけです。続いて、ここにデータを入れてみます。

lam = 1.0
sigma = 1.0
X = data[["x", "y"]]
y = data["c"]
label = data["c"] * 2 - 1
svm = NonLinearSoftMarginSVM(X , y, lam, sigma)
[w, b, Sup_index, alpha] = svm.fit()

ここまでは実はさほど計算時間はかかりません。決定境界をプロットしてみましょう。

def f(x, alpha, Sup_index, b, X, label):
  result = b
  for i in range(len(Sup_index)):
    j = Sup_index[i]
    result = result + alpha[j] * label[j]  * svm.kernel(X.iloc[j], x)
  return result

#色の設定
col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")

#メッシュグリッドを作成して決定境界をプロット
seq = np.linspace(0, 10, 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
for i in range(X1.size):
  Z.append(f(np.array([x1[i][0], x2[i][0]]), alpha, Sup_index, b, X, label))

plt.title("Non Lnear Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.0])
plt.scatter(data["x"], data["y"], c = col)

f:id:mt19:20181022155837p:plain
このように非線形の決定境界となります。SVM自体の話はここまでで、あとはどうやってモデルを良いものにしていくかが難しいところです。例えば今回はソフトマージンなので正則化係数、さらにはガウスカーネルなので分散\sigmaを事前に設定する必要があります。他には多項式カーネルであれば2つのハイパーパラメータを設定しなければなりません。クロスバリデーションなどを適用しながらこれらを適切に設定することが求められます。

参考

非線形の決定境界を書く際に参考になりました。
python で等高線を描くなら meshgrid して contour | コード7区

ソフトマージンSVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

前回はハードマージンSVMを実装しました。今回はそれの拡張のようなもので、ソフトマージンSVMを行ないます。ハードマージンSVMではデータは完全に線形分離できるものにしてくださいと書きました。
ハードマージンSVMによる分類 - lazy diary
しかし現実問題そんな都合の良いデータばかりではないでしょう。なんらかの理由でラベルの異なるグループに入ってしまうデータもあるはずです。ソフトマージンSVMでは、そのようなデータが入ることを許すよう条件を緩和します。しかし、どの程度まで許すかということが問題になります。そこで条件を緩和しすぎず、なおかつマージンを大きくすることを目指したのがソフトマージンSVMです。
実際に解析するにあたってその程度を決める指標として正則化係数を私たちは与える必要があります。これについてはクロスバリデーションで決める方法がおそらく最もシンプルだと思われます。

コード自体はハードマージンSVMとほとんど変わりません。二次計画問題で、条件が増える点と、サポートベクトル集合にも条件がさらに加わることだけ注意します。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class SoftMarginSVM():
  
  def __init__(self, X, y, lam):
    self.X = X
    self.y = y
    self.lam = lam
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

データを作ります。ハードマージンSVMと違い多少は入り乱れたデータでも平気です。

#create data
x_1 = 4.2
y_1 = 7
x_2 = 7
y_2 = 4.2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

これをSVMにかけてみましょう。

lam = 10.0
svm = SoftMarginSVM(data[["x", "y"]], data["c"], lam)
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.title("Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022142912p:plain
このように、ハードマージンSVMの厳しすぎる条件を緩和することでソフトマージンSVMでは分類を達成します。

ハードマージンSVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回はハードマージンSVMというものを実装します。
ソフトでもない上に非線形でもないという最も基本的なSVMです。ハードマージンSVMさえ理解してしまえば他のSVMはその拡張みたいなものです。

必要なものをimportしておきます。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

このcvxoptというのは二次計画問題のソルバーです。詳しくはこの記事の最後にあるブログを参考にしてください。google colabにもデフォルトではインストールされていないようなので、

pip install cvxopt

で導入しておきましょう。

データは相変わらずいい加減に作ったものです。とは言いつつも今回はハードマージンSVMですので、線形分離できないデータは困ります。かならず超平面で分離できるデータをチョイスしてください。2次元なら直線で完全に分離できるものです。

#create data
x_1 = 3
y_1 = 7
x_2 = 7
y_2 = 2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

ここからが本編です。流れとしては
1. データから二次計画問題を解く
2. パラメータwを求める
3. サポートベクトル集合を求める
4. サポートベクトル集合から切片bを求める
となっています。

class HardMarginSVM():
  
  def __init__(self, X, y):
    self.X = X
    self.y = y
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    G = matrix((-1) * np.identity(num_of_data))
    h = matrix(np.zeros((num_of_data, 1)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon:
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

これを上のデータに当てはめてみましょう。

svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022105408p:plain
黄色になっている部分はサポートベクトルです。SVMの特徴として実質的にサポートベクトルにしか依存しないので、次のような外れ値が入っても影響がありません。

data.loc[100] = [0, 20, 1]
svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022110158p:plain

参考

今回、二次計画問題を解くにあたって参考にしました。SMOアルゴリズムについてはまだ理解できていないため、今回はcvxoptというソルバーに頼りました。
二次計画法(Quadratic Programming)の概要とPythonソルバーcvxoptの使い方 - MyEnigma

XORのようなデータを非線形ロジスティック回帰で分類する

参考書は相変わらずこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回は、基底関数展開を利用した非線形のロジスティック回帰を行ないます。基底関数を求めて、データを変換、最後にロジスティック回帰を行なうだけです。今回は基底関数としてガウスカーネルを使いました。
非線形の恩恵を確かめるために、通常の線形ロジスティック回帰では分類できないXORのようなデータを今回は使います。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.cluster import KMeans
from sklearn import linear_model
from scipy.stats import multivariate_normal

cov = 0.5 * np.identity(2)
data1 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, -2.0), cov = cov, size = 25))
data2 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, -2.0), cov = cov, size = 25))
data3 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, 2.0), cov = cov, size = 25))
data4 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, 2.0), cov = cov, size = 25))

data1["c"] = 0
data2["c"] = 1
data3["c"] = 0
data4["c"] = 1

data = pd.concat([data1, data2, data3, data4], join = "inner")
data = data.reset_index(drop = True)
data = data.rename(columns = {0:"x", 1:"y", "c":"c"})

x = data["x"]
y = data["y"]
clusters = data["c"]

plt.scatter(x, y, c = clusters)

f:id:mt19:20181109161138p:plain
色付けをまともにしてないので見づらいですが、たしかにXORです。ここから前回と同じようにk-meansからの基底関数の構成およびデータの変換を行ないます。

#k-meansなどでクラスタリングする
dist = []
for i in range(1,30):
  km = KMeans(n_clusters = i)
  km.fit(data[["x", "y"]])
  dist.append(km.inertia_)
  
plt.plot(range(1,30), dist, marker = ".")

f:id:mt19:20181109161217p:plain
エルボー則に従って4くらいかなと考えられます。

from sklearn import linear_model
#クラスタ数を4とする

#ガウス型基底関数
def gauss_func(x, mu, sigma):
  return np.exp((np.sum(np.power(x - mu, 2))) / (-2 * sigma))

m = 4
km = KMeans(n_clusters = m)
km.fit(data[["x", "y"]])
cluster_center = km.cluster_centers_
dist = np.power(km.transform(data[["x", "y"]]), 2)
cluster_num = km.predict(data[["x", "y"]])
least_dist = np.zeros(shape = (m, 1))
for i in range(len(data)):
  for j in range(m):
    if cluster_num[i] == j:
      least_dist[j] += dist[i][j]
for i in range(m):
  least_dist[i] = least_dist[i] / np.sum(cluster_num == i)

transformed_X = pd.DataFrame(np.zeros(shape = (len(data), m)))
for i in range(len(data)):
  for j in range(m):
    transformed_X.iloc[i, j] = gauss_func(data[["x", "y"]].iloc[i,:], cluster_center[j], least_dist[j])

さて、カーネルパラメータの推定とカーネルを使ってデータの変換まで終わったのでロジスティック回帰を実行します。

#ロジスティック回帰の実行
log_reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga", C = 1000.0, max_iter = 1000).fit(transformed_X, clusters)
print(log_reg.score(transformed_X, clusters)) #->1.0

とりあえず分類できました。決定境界をプロットしてみます。

seq = np.linspace(data["x"].min(), data["x"].max(), 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
testdata = np.hstack((x1, x2))
transformed_testdata = pd.DataFrame(np.zeros(shape = (len(x1), m)))
for i in range(len(x1)):
  for j in range(m):
    transformed_testdata.iloc[i, j] = gauss_func(testdata[i,:], cluster_center[j], least_dist[j])
Z = log_reg.predict_proba(transformed_testdata)

plt.title("Non Linear Logistic Regression")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z[:,0]).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.5])
plt.scatter(data["x"], data["y"], c = clusters)

f:id:mt19:20181109161626p:plain
ガウスカーネルっぽい丸みを帯びた決定境界です。
基底関数になるものならばどれを使ったって構いません。(うまく分類できるかは別ですが)。例えばx^n(n = 0, 1, 2,\dots)を用いても良いです。

基底関数展開を使った非線形回帰

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

この第3章非線形回帰モデルの内容です。

今回は基底関数を使った非線形回帰です。
ここで言う基底というのは線形代数で学ぶ基底と同じ意味合いです。ここでは基底"関数"ですので、その基底の要素が関数になっているわけです。

非線形回帰モデルで最も有名な例としてd次元多項式回帰を考えると、その基底関数は\{1,x,x^2,x^3,\dots,x^d\}となります。この関数にデータx_iを入力してd次元のデータを作り線形重回帰モデルを学習させることで多項式回帰は非線形回帰として機能しています。

作用させる関数を変えればより表現能力の高いモデルが可能になるのではということで基底関数として他にも様々なものが考えられています。最も典型なのがガウス型基底関数で

\displaystyle{
\phi_j(x) = \exp{\left (-\frac{\|x - \mu_j\|^2}{2h_j^2}\right )}
}
という関数になります。パラメータとして位置を決める\muとばらつきを決めるhがあります。これはデータから推定します。今回はこの基底関数を用いることにします。具体的には次の手順でモデルを作ります。

1. 基底関数の個数を決めるためにデータに対してk-meansなどでクラスタリングを行なう。このときのクラス数をnとしておく。
2 各クラスに属するデータの平均と距離による分散を求め、合計nクラスぶんのガウス型基底関数を構成する。
3. データにn種類のガウス型基底関数をそれぞれ作用させ、そのデータに対して線形重回帰モデルを学習させる。

この手法は多次元にも適用できるので、今回は可視化することを考えて説明変数は2次元にしておきました。

from sklearn import linear_model
from sklearn.cluster import KMeans
from scipy import stats
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def create_multidata():
    N = 20
    residual = stats.norm.rvs(N) / 10
    x = np.linspace(-5, 5, N) + residual
    y = np.linspace(-5, 5, N) + residual * 2.0
    return x, y, (1 / 10) * x**2 - x + y + np.sin(x * y * 10) + residual


def gaussian_basis(x, mean, std):
    result = np.empty((x.shape[0], 1))
    for i in range(x.shape[0]):
        temp = (x[i, :] - mean).reshape(-1, 1)
        result[i] = np.exp(-np.linalg.norm(temp) / (2 * std**2))
    return result.flatten()


def gaussian_basis_func_reg(x, y, n):
    kmeans = KMeans(n_clusters=n).fit(x)
    label = kmeans.predict(x)
    mean = np.empty((n, x.shape[1]))
    std = np.empty((n, 1))
    transformed_x = np.empty((x.shape[0], n))
    for k in range(n):
        x_k = x[label == k]
        mean[k, :] = np.mean(x_k)
        std[k, 0] = np.sum([np.linalg.norm(x_k[i, :] - mean[k, :]) for i in range(x_k.shape[0])]) / x_k.shape[0]
        transformed_x[:, k] = gaussian_basis(x, mean[k, :], std[k, 0])
    model = linear_model.LinearRegression().fit(transformed_x, y)
    return model.predict(transformed_x)


if __name__ == "__main__":
    x, y, z = create_multidata()
    X = np.vstack((x, y)).T
    Y = z.reshape(-1, 1)
    model = linear_model.LinearRegression().fit(X, Y)
    non_linear_pred = gaussian_basis_func_reg(X, Y, n=3)
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    ax.scatter(x, y, z, label="observed", c="blue")
    ax.plot(x, y, model.predict(X).flatten(), label="linear", c="red")
    ax.plot(x, y, non_linear_pred.flatten(), label="non linear", c="green")
    ax.grid(True)
    ax.legend(loc="upper left")
    plt.show()

f:id:mt19:20190817093312p:plain

モデル選択はモデルの有効自由度を求めてAICを使うか、クロスバリデーションがあるようですが後者のほうが実装しやすいかなとは思います。

ガウス型基底関数以外にはthin plate spline関数とか逆多項式型関数とか色々あるようですが、どれもこれも自分で実装しなければならない(sklearnに無い)上にパラメータの推定法についてもちょっと調べないといけない状況です。あまり人気がないんですかね(カーネル法を使った方が手っ取り早い?)。

今回は以上です。


更新履歴
2018/10/11 記事作成
2019/08/17 復習するついでに1次元だった説明変数を多次元に拡張しました。