バナナでもわかる話

開設当初は計量経済学・統計学が専門の大学院生でした。今はデータを扱うお仕事をしています。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

機械学習のバイアスバリアンス分解を古典的な統計学から丁寧に解説する

今日は、バイアスバリアンス分解の話をしようと思います。


目次


スポンサーリンク



はじめに

「バイアスバリアンス分解」と言う言葉を聞いて、皆さんは何を思い浮かべるでしょうか。

恐らく古典的な統計学だけを勉強した人であれば、パラメータの推定量に関するMSEの分解を思い浮かべるかもしれません。

一方機械学習から勉強した人であれば、トレーニングデータで評価した予測値と予測変数(目的変数)とのMSEを思い浮かべることでしょう。

実は、後者のことを一般にバイアスバリアンス分解と呼び、前者はMSEをバイアスとバリアンスに分解できるだけの話で特に名前は与えられていない(と思われる)のですが、当然、両方MSEですし、やってることは同じなので、両方勉強している人であれば、そんなもの分けて議論することがナンセンスだともなることでしょう。


しかし、恐らく機械学習に馴染みのない方からすると、機械学習サイドでよく言われるところのバイアスバリアンス分解はわかりにくく、

躓きポイントじゃないかな〜と私は考えています。



そこで今回の記事では、まずはじめに古典的な統計学で語られるところのMSEの分解とその理解について説明した上で、それとのモチベーションの違いに触れながら機械学習で語られるところのバイアスバリアンス分解について説明しようと思います。

そして、最後にオマケとして、最近話題のdouble descentの話にも触れようと思います。

古典的な統計学におけるMSEの分解

パラメータの推定

まずはじめに、古典的な統計学の目的を確認しておきましょう。
統計学では、あるデータXの性質を'分布'という観点から考察したいというモチベーションから出発します。
ここであるデータXの性質とは「期待値や分散の値」だけでなく、「Xに影響を与える変数は何か」なども含みます。

とにかく、データXの分布(当然この分布は仮定する)のパラメータを推定し、それを解釈利用することでXの性質を述べたり、その述べた性質が妥当か検証したり、パラメータの値を使って予測したりといったことを行いたいわけです。


なので、とにもかくにも仮定した分布のパラメータの値を、持っているデータを使ってうまく当てたいわけです。


ここでXに対して、何かしらのパラメータ \thetaによって決まる分布を仮定します。

パラメータ \thetaの値は神様しかわからないわけですが、持っているデータXを用いて得られる推定量 g(X)で何とか誤差が小さくなるようにパラメータの値を推定したい!と言う話になってきます。


そこで、次のような評価方法Lを考えることになります。


 L=(\theta-g(X))^2

このLは神様しか知らない真の値と推定に用いる値との二乗誤差です。

Xは確率変数を仮定しており、確率的にばらつくはずなので、このLをXの期待値で評価すると

 MSE=E[(\theta-g(X))^2]


これが、MSE (Mean Squared Error) です。

と言うわけで古典的な統計学でMSEを用いる基本的なモチベーションとしては、


あるデータXの性質をパラメータに置き換えて述べたい
しかし、パラメータの値はわからないのでうまい推定量で推定したい
うまい推定量とは何かを検討するにあたり、MSEで評価したい


ということになります。

MSEの分解

さて、MSEですが、推定量のバイアス分散に分解することができます。

ちなみにここでいう推定量のバイアスとは何を指すかというと、
推定したい真の値と、推定量の期待値の差です。
つまり \theta-E[g(X)]がバイアスです。


それを示していきましょう。

 MSE

 =E[(\theta-g(X))^2]

 =E[(\theta-E[g(X)]+E[g(X)]-g(X))^2]

 =E[(\theta-E[g(X)])^2+(g(X)-E[g(X)])^2+2(\theta-E[g(X)])(E[g(X)]-g(X))]

ここで、最後の項は期待値をとるとゼロになりますね。
 E[2(\theta-E[g(X)])(E[g(X)]-g(X))]=0
Xの期待値を考えていることに注目してください。


よって、
 MSE

 =E[(\theta-E[g(X)])^2]+E[(g(X)-E[g(X)])^2]

 =(\theta-E[g(X)])^2+V[g(X)]

1項目はバイアスの二乗ですね。
2項目は推定量の分散です。

こうしてMSEをバイアスと分散に分解することができました。


ここまで見れば、例えば不偏推定量(バイアスが0の推定量)を考えるモチベーションもわかるかと思います。
MSEを小さくしたいが、MSEは正確には分散とバイアスの二つの要因から成る。
バイアスが0の推定量のクラスであればとにかく分散を小さくすることだけに注力すればいいので、考えやすいですね。

当然、モチベーション次第では、若干バイアスが生じてもそれ以上に分散を小さくできるような推定量があれば、不偏である必要もないですね。
そんなことを考えることができるわけです。


機械学習におけるバイアスバリアンス分解

予測と分類

一般に機械学習と言われる手法では、上記の話とはモチベーションが少しだけ違います。
パラメータには重きを置きません。


とにかく目的変数を予測(分類)したい!!とにかく正確に予測したい!!というモチベーションで話が進みます。

いやいや、先ほどの統計学の説明でも予測はモチベーションの一つとしてあったじゃないかというのは当然の指摘ですが、先ほどの古典的な統計学のモチベーションのメインはとにかくパラメータにあります。その成果を使って予測を行うこともあるというだけです。


古典的な統計学の文脈で良しとされる手法は、予測変数のMSEを出来うる限り小さくしているわけではないのです。

一方、機械学習でよしとされる手法は、パラメータのMSEを最小化しているわけではないのです。


若干ミスリードの感があるので補足しておくと、これは、統計学は予測に興味がないといっている訳ではありません。統計学は予測に"も"興味はあるが、他のことにも興味があるため、パラメータ目線で話を進めることが多かったというだけです。そのため機械学習を予測に特化した統計学として解釈することも当然可能です。




では、機械学習の話に戻ります。
例えばデータXを使ってデータYの値を予測したいといった状況を考えることにしましょう。

この時、得られたデータ X^{(T)}を使って g(X)というYを予測するための関数(予測器)を考えて、新しいXが与えられた時の'未知のY'を予測したいというのが機械学習のモチベーションです。

ここで、添字 (T)を与えてモデルを作る時に使っている X^{(T)}と、予測したいYに対応するXを区別していることに注意してください。

既に持っているデータは、答えであるYがもう分かっているのだから、わざわざ予測する必要なんてありません。
まだわからないものを予測してこその予測器です。

オーバーフィット

既に持っているデータ X^{(T)}とそれに対応する答え Y^{(T)}の組をテストデータと呼びます。

実はテストデータの二乗誤差

 \sum_{i=1}^{N}(Y_i^{(T)}-g(X_i^{(T)}))^2

を0にする方法は簡単で、パラメータの次元を増やす(モデルを複雑にする)ことで誤差0の理想のモデル(!?)が完成します。


しかし、これには問題があって

このようにテストデータに対する誤差最小化モデルは、新しいデータに対する予測性能がありません

これは当然で、テストデータの二乗誤差を最小化するモデルを作っただけで、新しいデータに対する二乗誤差最小化モデルを考えていないからです。


このようなテストデータにだけ過剰にフィットしている状況をオーバーフィットと呼んだりします。

よくこんな図とともに説明されます。
f:id:bananarian:20190917234218j:plain*1
点線がテストデータで、横軸のcomplexityは複雑度ですが、とりあえずパラメータの次元だと思ってください。縦軸もRiskですが、今回の話で言うところの二乗誤差だと思ってもらえれば問題ないです。

実線が新しいデータで見たものです。複雑度が増すと一定のラインを超えると、Riskが猛烈な勢いで増えています。


じゃあどうすればオーバーフィッティングを回避できるのかというと、評価の基準をテストデータに関する二乗誤差ではなく、新しいデータに関する二乗誤差にしてやれば良い訳です。


バイアスバリアンス分解

それでは、MSEを考えましょう。改めて文字を整理します。


予測器を作るにあたって用いるテストデータセットをDとして次のように書くことにします。

 D= \{(X_1^{(T)},Y_1^{(T)}),...,(X_N^{(T)},Y_N^{(T)})\}


このテストデータセットDを用いて作った予測器gを次のように書くことにします。

 g(x;X^{(T)},Y^{(T)})=g(x;D)


これは、gは新しい予測のための変数xを入れることでxに対応する予測値を出力する関数であるが、
作るにあたって得られたテストデータDの値を用いている(つまり、Dを条件付けたもとで得ている)ことを明示的に示すために、このような書き方をしています。


さて、テストデータとは区別して、データXが与えられた時に生じるであろう値Yをg(x)を使って予測したいので、次のような誤差Lを持って評価すれば良いことがわかりますね。


L=(Y-g(X;D))^2


ただしここで、Yは確率変数です。そのため、期待値をとって次のように考えてやる必要があります。

 E_{y|x}[(Y-g(X;D))^2]

演算子 E_{y|x}はxを条件づけた時のyの期待値を求めていることを表します。ちゃんと積分の形に直すと次のようなことを言っているのと同じです。

 \int (Y-g(X;D))^2 p(y|x)dy


さて、実はもう一つ確率変数があります。Yを予測するにあたって、Xは固定なのですが、データセットDは固定していません

モデルを考えるにあたって確率変数の組であるトレーニングセットを使って予測器を作っているため、この g(X;D)は確率変数です。
要はXを固定しているだけで X^{(T)}を固定した訳じゃないのです。よって、こいつの期待値も取る必要がありますね。

 MSE=E_{D|x}[E_{y|x}[(Y-g(X;D))^2]]


これが機械学習で考えているところのMSEです。

式は一緒でも、考えているモチベーションが若干異なるため、このような形になる訳ですね。
さて、分解してみましょう。手順は先ほどと同様です。


 MSE

 =E_{D|x}[E_{y|x}[(Y-g(X;D))^2]]

 =E_{D|x}[ E_{y|x}[ (Y-E_{y|x}[Y]+E_{y|x}[Y]-g(X;D) )^2 ] ]

 =E_{D|x}[ E_{y|x}[ (Y-E_{y|x}[Y])^2 + (g(X;D)-E_{y|x}[Y])^2 ] ]

若干複雑ですが、中身をしっかりみてください。右の項はYの値によらずバイアスの二乗の形であることがわかりますね。


さて、続きの計算を行いましょう。
 =E_{D|x}[ V_{y|x}[Y]+(g(X;D)-E_{y|x}[Y])^2]

 = V_{y|x}[Y]+E_{D|x}[(g(X;D)-E_{y|x}[Y])^2]

さて、右側の項をさらに分解しましょう。分解方法は今までの手順と同じです。


  =V_{y|x}[Y]+E_{D|x}[(g(X;D)-E_{y|x}[Y])^2]

 = V_{y|x}[Y] + E_{D|x}[(g(X;D)-E_{D|x}[g(X;D)]+E_{D|x}[g(X;D)] -E_{y|x}[Y])^2]

 =V_{y|x}[Y] + E_{D|x}[(g(X;D)-E_{D|x}[g(X;D)])^2+(E_{D|x}[g(X;D)]-E_{y|x}[Y])^2]

 =V_{y|x}[Y] +V_{D|x}[g(X;D)]+(E_{D|x}[g(X;D)]-E_{y|x}[Y])^2

機械学習モデルにおいては、しばしばこの式の第三項目をバイアスと呼んだりします。
第二項目はバリアンス、そして第一項目はトレーニングデータセットと無関係なデータの分散です。

この第一項目についてはトレーニングデータではどうしようもないので、いじれません

よって、予測においてはトレーニングデータとモデルを使っていじることが可能な第二、第三項が重要になります。


まとめ

結局、やっていることも、評価方法も同じであるということが分かっていただけたと思います。

分析において評価したいものが違っていて、それは手法の背後にあるモチベーションからこそなのです。

そして、評価したいものが違うがゆえに、期待値演算が二重になって若干複雑になっているのが機械学習(or 最新の統計学)におけるバイアスバリアンス分解なのです。


double descent

先ほど、オーバーフィッティングで次のような図を確認しました。
f:id:bananarian:20190917234218j:plain

しかし、2018年、実は次のような状況にあることが[Belkin, Mikhail, et al., 2018]によって確認されました。
f:id:bananarian:20190917234853j:plain*2

実は、もっとモデルの複雑度を増やすとまたRiskが減少し始めるらしい。これをdouble descentと呼んだりします。
これは、古典的な統計学の手法では、そもそもこの先がわかると言う状況がなかったので知られていなかったのですが、例えばニューラルネットワークなどを考えるとこのような状況が発生します。

この仕組みはまだ完全に判明していないので、理屈を説明することが出来ないのですが、

これは今注目を浴びているトピックで、テストデータでRMSEの谷を見るだけではどうやらダメらしいということで研究が進んでいる最中だったりします。


*1:Belkin, Mikhail, et al. "Reconciling modern machine learning and the bias-variance trade-off." arXiv preprint arXiv:1812.11118 (2018)より

*2:Belkin, Mikhail, et al. "Reconciling modern machine learning and the bias-variance trade-off." arXiv preprint arXiv:1812.11118 (2018)より

単純パーセプトロンとマージンに関するまとめ

今日は単純パーセプトロンの話をします。
実際の現場で役に立っているのは多層パーセプトロンや、SVMですが、その基礎となるパーセプトロンについて理解しておくことは重要ですので、とりあえず一通りまとめていきます。


目次

スポンサーリンク


単純パーセプトロンとは何か


まず、単純パーセプトロンとはなんですか?って話をしておきます。

簡単にいうと、「直線(超平面)を引っ張って、空間を二つに分けることで、2クラス分類をやろう」というのがこの手法にあたります。
簡単な例を示すとこんな感じ。

f:id:bananarian:20190802215915p:plain

2つの領域に分けて、うまいこと2種類のクラスに分類しています。

直線を引っ張るだけなら、線形回帰か??と思われる方もいらっしゃるかもしれませんが、線形回帰とは違います。その辺の話も後でしようと思います。


数学的な理解


まず、適当な直線を考えましょう。未知パラメータベクトル \omega=(\omega_0,...\omega_k)、特徴量(統計の文脈で言うところの説明変数) X_1,...X_kを用いて次のような直線を考えます。

 y=\omega_0+\omega_1 X_1+\cdots+\omega_k X_k

気持ちとしては、直線の上側と下側でクラスが分かれると嬉しいので、次のような関数 \phiをに突っ込みましょう。


 y≥0の時
 \phi(y) = 1

 y<0の時
 \phi(y) = -1


これで、yが実数範囲でどのような値を取ろうとも、1か-1かどちらかに分類することができます。
ちなみに、ここで1と−1にしている理由もしっかりあります。これは後々の解説で説明します。




さて、 ここで y≥0ですが、

 y=\omega_0+\omega_1 X_1+\cdots+\omega_k X_k

でしたので、

 \omega_0+\omega_1 X_1+\cdots+\omega_k X_k≥0、つまり \omega_1 X_1+\cdots+\omega_k X_k≥-\omega_0

変形後の左側の式は、切片0の直線です。このように見方を変えると、 \omega_0はちょっと他の係数とは役割が異なり、領域の閾値を表していることがわかります。
そのため、機械学習の文脈では、この切片のことを特に他の係数と区別して、バイアス項と呼んだりもしています。


ここまで読めば、線形回帰との違いはわかりますよね。

数式的な観点で言えば
そもそも線形回帰とは違って誤差項を仮定しているわけではありませんし、確率的な変動も考えていません。
あくまで、2クラスに分類できると仮定して、直線を引っ張っているだけであることに注意してください。
(当然確率的な変動を考えて定式化することもできる)

モチベーション的な観点から見ても
線形回帰はデータ点が確率的な変動を伴いながら線形に並ぶことを想定した統計手法ですが、

単純パーセプトロンはデータを特徴量空間に並べた場合に、直線で二分することができるということを想定した分類アルゴリズムです。


パラメータの探索

パラメータ探索イントロ

さて、整理したところで、次にパラメータの値をどう決めようと言う問題にシフトすることにします。


統計学にせよ、機械学習にせよ、どんな値が尤もらしいかを考えて、パラメータの値を決めなければなりません。

言い方を変えると、パラメータの値を見つけるためには、妥当な計算しやすい評価方法を考え、その方法を持って評価する必要があります。


さて、ここで次のようなデータが与えられているとしましょう。

f:id:bananarian:20190802225239p:plain


正解ラベルとしてLabel(1 or -1)が与えられているデータです。このデータを使って、2クラス分類のアルゴリズムを単純パーセプトロンを作り、例えば新しいデータがやってきたときに、特徴量を使って正しく分類したいと言ったことがモチベーションです。


そこで、まず問題になるのは


どんな値で評価しよう??と言うことです。

パラメタ探索の評価

すぐ思いつくのは、

予測ラベル \hat{\phi(y_i)}と、正解ラベル  Label_iを用いて評価関数Eを次のように設定すれば良いのではないか?ということです。

 E = \sum_{i=1}^n (\hat{\phi(y_i)} - Label_i)^2

こうすれば、 \hat{\phi(y_i)} Label_iが一致すれば0,不一致ならば 2^2となり、全て正しく分類できている場合に限り E=0となります。

Eが極力小さくなるような \omegaを探せばいいじゃないか!めでたしめでたし


とはなりません。このEには欠陥があります。
Eは離散値しかとらないので、微分を使った最適化法が行えません。
とすると適当な離散最適化から攻めたり、データが少なければ総当たりという方法も考えられますが、やはり離散値では使い勝手がよくありません。


そこで、うまく離散値をとるような評価関数を考えます。


次のようなEを考えてみましょう。

 E = -\sum_{i=1}^n y_i Label_i 1_{(y_i Label_i<0)}

ただし
 y_i = \omega_0+\omega_1 X_{(1,i)}+\cdots+\omega_k X_{(k,i)}
 1_{(y_i Label_i<0)}は定義関数(指示関数)で、 y_i Label_i<0の時に1、そうでないときに0をとる関数です。

ちなみに、 y_i Label_iは各点における関数マージンとも呼ばれます。



こうすることで、誤差が出る場合には、Eは値が増え、正しい場合には0のまま、しかも連続的な値をとる関数と見ることができます。

このEを最小化(0に近づける)することを考えます。


パラメータの探索(確率的勾配法)

では、このEを \omegaで微分してみましょう。

 \frac{\partial E}{\partial \omega_s} = -\sum_{i=1}^n X_{(s,i)} Label_i 1_{(y_i Label_i<0)}

ここで、勾配が0であれば適当な局所解であるはずです。

よって、とりあえず適当に \omegaの初期値として \omega_0^{(0)},...\omega_k^{(0)}を与えて、そこから勾配に従って修正を行なっていきます。

つまり

 \omega^{(t)} = \omega^{(t-1)} + \frac{\partial E}{\partial \omega}

これを繰り返していくことで、適切に分類できるパラメータ値に近づけていきます。


ここで、 Label_i 1_{(y_i Label_i<0)}は次のように書いている本もあったりしますが、意味としては同じですね。

 (y_i - Label_i)

こっちの方が見やすいので、以降は次のように書くことにします。

 \frac{\partial E}{\partial \omega} = -\sum_{i=1}^n X_{(s,i)} (y_i - Label_i)= \sum_{i=1}^n X_{(s,i)} ( Label_i-y_i)


ということで、これで、探索アルゴリズムも完成しました。めでたしめでたし。


より効率の良い又は精度の良いアルゴリズム

例えば、 X_1が20~30の範囲をとるような場合を考えます。この時、 \omegaの更新は20ほどの値を足したり引いたりして近づいていくことになりますが、この更新の値幅って20より2の方が正確だし、2より0.2の方がより正確な値を探すことができるはずですよね。

ただし一方で、この値が細かくなればなるほど収束時間が長くなるはずですね。

つまり、この更新幅は調整出来たようが嬉しいわけです。


ということで、先ほどのアルゴリズムを少し改善して、次のように学習率 \eta \in [0,1]を与えてみましょう。


 \omega^{(t)} = \omega^{(t-1)} + \eta \frac{\partial E}{\partial \omega}


この学習率を適当な値に変えて、収束するパラメータを探していくというアルゴリズムになるわけです。


スポンサーリンク



マージンに関する話

関数マージン

データ点関数マージン

先ほど出てきたマージンについてもう少し考えてみることにします。

今、与えられているデータとして (X_i,Label_i)がありますが、
※特徴量は以下 X_iというベクトルの形でまとめて表記し直す。
 X_i = (X_{(1,i)},...X_{(k,i)})^{'}

この時、あるパラメータ \omega=(\omega_1,...\omega_k)とバイアスパラメータ \omega_0を使って、次のような直線(超平面)を引くのでした。


 \omega X_i +\omega_0

さらに、先ほどちらっと出てきましたが次のような \gamma_iデータiに対する関数マージンと呼んでいました。

 \gamma_i = Label_i (\omega X_i +\omega_0)

超平面関数マージン

さて、仮にパラメータ \omega,\omega_0を適当な値に固定したとします。
そうすると、超平面は適当な値に固定されるわけですが、この場合、データ点で見たときの関数マージンはデータ点の個数分(つまり、今回ならn個)出てきます。

 \gamma_1,\cdots,\gamma_n

このうち最小の値をとる関数マージン min\{\gamma\}を超平面における関数マージンと呼びます。

幾何マージン

データ点幾何マージン

関数マージンだと、正直数式的には、そうですかと言った感じですが、直観的に何をやっているのかがわかりにくい。
そこで、次のようにデータ点における関数マージンを正規化してやることにします。


 \frac{\gamma_i}{||\omega||}

これを、幾何マージンと呼びます。
分母はバイアスパラメータ以外のパラメータに関するノルムです。

これのどこがわかりやすいのかというと、絶対値を取れば、点と直線の公式になっていますよね?

つまり、正解率100%の超平面を引いた場合、幾何マージンは超平面と点の距離になっているということになります。
心なしかイメージしやすくなりました。

超平面の幾何マージン

関数マージンの場合と同様ですが、各データ点における幾何マージンのうち、最小の値をとる幾何マージンを、超平面における幾何マージンと呼びます。

最大マージン超平面

ここまでは、データ点を固定して、各点におけるマージンを見たり、超平面を固定したときに、特に小さいマージンを見たりしましたが、

実際の超平面は様々考えられます。


考えられうるすべての超平面における幾何マージンを考えてみましょう。
当然、分類に誤りがある超平面であれば、幾何マージンは負の値を取ります(各データ点における最小の幾何マージンが超平面の幾何マージンでした)。

一方、100%正解する超平面の複数考えられますが、その場合、一番近いデータ点からの距離が近ければ幾何マージンは小さくなりますし、離れていれば幾何マージンは大きくなります(実際にフリーハンドで描いてみるとわかりやすい)。


100%正解するのであれば、データ点スレスレの直線を引くよりも、余裕を持った直線を引いた方が良いですよね。

つまり、超平面の幾何マージンは大きい方が嬉しいわけです。

そこで、すべての、超平面の幾何マージンのうち、最大の値をとる幾何マージンを最大マージン超平面と名付けます。

当然ですが、線形分離可能なデータであれば、この最大マージン超平面は正の値を取ります。

重要定理について

さて、以下は一通り重要な定理を確認していきます。
共立出版のサポートベクターマシン入門が非常によくまとまっていて、そのまま証明を書くと丸写しになってしまうので、証明自体は本に譲ります。
ここでは、上の説明と結びつけながら、定理の意味について説明していきます。

若干和訳に難がありますが、非常に良い本ですので、副読本にオススメです。

Novikovの定理

ある、2クラスのデータが両方含まれるデータ集合(トレーニング集合)に対して、線形分離可能を仮定する。
更に次の値 Rを考えます。

 R = max_{1≤i≤n} ||X_i||

ここで、すべてのiについて、 ||\omega_{opt}||=1で、更に

 Label_i (\omega_{opt} X_i+\omega_0)≥\gamma

※ただし \gammaは既に与えている分離超平面におけるマージン

が成り立つような \omega_{opt}が存在すると仮定します。
つまり、より良い超平面(または同等)が存在するという仮定です。


このとき、このデータ集合について、パーセプトロンアルゴリズムが失敗する数は高々

 (\frac{2R}{\gamma})^2

に抑えられることがわかります。


これがNovikovの定理です。


気持ちとしては、線形分離可能であるならば、アルゴリズムが失敗する回数の上限がこの数で抑えられるので、高々有限回の試行で収束するという保証です。

更に、気になる点としては、失敗の上限は学習率に依存しないというところですね。


Freund and Schapire の定理

マージンスラック変数

この定理では、線形分離不可能な場合はどう考えるんだという疑問に対する解答が与えられます。
線形分離不可能な場合、マージンで考えるのは難しいので、新たな概念としてマージンスラック変数 \xiを導入します。

超平面に対するマージンを \gamma>0として固定する。これを以下ターゲットマージンと呼ぶ。
このとき、各点が超平面からこのターゲットマージンを持つにあたって、どの程度失敗しているかを \xi_iとして定量化する。

 \xi_i((X_i,Label_i), (\omega,\omega_0),\gamma) = max\{0, \gamma-Label_i (\omega X_i +\omega_0) \}

これをマージンスラック変数と呼ぶことにします。

定理詳細

重複もなく、2クラスのデータが含まれたデータ集合に対して、||X_i||≤Rとする。
更に (\omega,\omega_0) ||\omega||=1を満たす任意の超平面とする。
 \gamma>0

このとき、変数Dを次のように定義する。

 D = \sqrt{\sum_{i=1}^n \xi((X_i,Label_i),(\omega,omega_0),\gamma )^2}


このとき、パーセプトロンアルゴリズムの最初の1回目の実行での失敗の数は高々

 (\frac{2(R+D)}{\gamma})^2

になる。



つまり、線形分離不可能なデータ集合では、最終的にアルゴリズムの結果が振動して、収束しないわけだが、そうは言っても最初の1回の試行における失敗の数は上限があると、この定理から述べることができるわけです。


今回の記事の話に関する機械学習関連のおすすめ図書

サポートベクターマシン入門

ITエンジニアのための機械学習理論入門

パターン認識と機械学習 上

Python機械学習プログラミング

統計検定1級対策解説~情報量基準編~

今回はAICです。
AIC周りの話は、カルバックライブラー情報量とAICにどういう結びつきがあるのか概要を知っておくことが重要になってきます。

スポンサーリンク



AICとは

AICは次のような式で表されます。


 AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)


まず \hat{\theta}は、n個のサンプルから得られるパラメータの最尤推定量です。そして、 dim(\theta)は分布の自由度となっています。

このAICとはなんぞやというのを理解することがこの章で重要な話になってきます。

何の役に立つんだって話

このAIC、何の役に立つのかというと、オーバーフィットを回避するためのモデル選択に利用できます。

オーバーフィットとは

例えば、こんな感じのデータがあったとします。
f:id:bananarian:20190515123427p:plain

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)

で、どういう関数を当てはめたら良いかよくわからないので、最小二乗法(又は最尤法)を使って、多項式を当てはめてみることにするとします。


多項式っていうのはこういうやつですね。

 f(x)=\sum_{m=1}^M \alpha_{m-1}X^{m-1}+\epsilon

これはM-1次の多項式です。
Mはどうやって決めるのか問題がありますが、最尤法であれば、対数尤度最大、最小二乗法であれば二乗誤差を最小にするように、Mを適当に動かしてやると探すことが出来ますね。ちなみに、機械学習界隈ではよくRMSEが利用されますがどれを使おうとさしたる違いはありません。


で、普通に最小二乗法、又は最尤法を使うと次のような曲線が得られます。
f:id:bananarian:20190515124216p:plain

「おぉ~~全ての点を通ってて良い感じ~!」

と満足してはいけません。この曲線には問題があります。

同じ分布から新しいサンプルをこんな風に取ってきたとします。

xtest=seq(0,1,0.01)
test1=sin(2*pi*xtest)+rnorm(length(xtest),sd=0.3)

今の曲線をこの新しいサンプルにあてはめるとこんな感じ。
f:id:bananarian:20190515124809p:plain

変なところにすっ飛んでいってる所があり、予測が見当違いになってしまいます。


このように、既存のデータに極端にフィットしてしまい、新しいデータに対する予測性能が無くなってしまう現象をオーバーフィッティングと呼びます。

パラメータの次元(今回ならM)も最適化しようとした場合に、最尤法や最小二乗法を用いると、このような状況が発生します。それを回避するために利用される手法がAICです。

AICの式の意味

AICの式をもう一度確認してみます。

 AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)


「このAICを最小化するようなパラメータを選択するのが良い」というのがAICが利用される場合の主張です。

第1項は最大対数尤度ですね。普通の最尤法ではこのAICの第一項を最小化(マイナスがついているので)するような値を選択しているわけですが、それではパラメータの次元を弄り始めるとさっきのような問題が発生します。

そこでAICでは第二項目としてパラメータの自由度を加えます。

つまり、パラメータの自由度を増やせば増やすほど罰則がかかるようになっているわけですね。

AICを使ってみる

では、自由度をいじりながらAICの値を確認してみます。簡単にプロットしてやると次のようになりました。

f:id:bananarian:20190515130147p:plain

最小点は6ですね。つまり、「10次多項式なんていらない、6次で十分だ」とAICで主張できるわけです。

ちなみに6次でさっきのプロットをしてみるとこんな感じ。
f:id:bananarian:20190515130703p:plain
f:id:bananarian:20190515130715p:plain

見た目的にも最小二乗法を使った場合よりも良いですね。

余談(機械学習でよく使われる手法の話)

ちなみに、機械学習や予測の文脈では次のような手法も使われます。

持っているデータを

・最小二乗法(最尤推定)用のトレーニングデータ
・次元Mを決めるためのテストデータ

にわけることで、先ほどの問題を解決するといった手法です。


こんな感じでテストデータとトレーニングデータがあったとします。

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)
test=sin(2*pi*x)+rnorm(11,sd=0.3)

この2つのデータに対して、トレーニングデータで学習した最尤推定量を元に、パラメータの自由度をいじって、その際のRMSEを見ることを考えます。

すると次のような推移になります。
f:id:bananarian:20190515131537p:plain

テストデータを用いた場合は、途中でRMSEの減少が止まっています。
つまり、テストデータを用いた場合、10次元も6次元もいらない。3次多項式で十分だと主張できることになります。

先ほどのAICは、尤度を利用していることからも分かる通り、サンプルサイズが少ない場合十分に機能しない可能性があります。そのため、データセットをこのように分ける手法が取られることもあるわけです。


カルバックライブラー情報量

カルバックライブラー情報量は、2つの分布間の(ある種の)距離を測る考え方です。
次のように定式化します。

2つの密度関数f,gについて

 KL(f|g)=\int f(x)log(\frac{f(x)}{g(x)})dx

2つの密度関数の対数の差の期待値を取っているのが見えますか?
もし、二つの密度が全く同じであれば0になります。

カルバックライブラー情報量とAICの関係性

ある分布 f(y;\theta)から標本Xが得られたとします。
この標本を使って分布のパラメータ \thetaを推定することを考える場合、その推定量 \hat{\theta}を使って尤度 f(y;\hat{\theta})を考えた場合、真の分布 f(y;\theta)に近くならなければならないはずです。


ここで、未知の標本(架空の物、得られていない) Yを使って真の分布と推定量を突っ込んだ分布の距離をカルバックライブラー情報量で考えてやると次のようになります。


 KL(f(y; \theta)|f(y;\hat{\theta}))=\int f(y;\theta)log(f(y; \theta))dy -\int f(y; \theta) log(f(y; \hat{\theta)}) dy

この値を小さくするためには、第1項目は弄りようが無いので、第2項目を考える必要が出てきますね。
つまり、第二項目の

 \int f(y; \theta) log(f(y; \hat{\theta)}) dy

を最大にするような推定量が汎化性能(予測とかの意味で良い性質を持つこと)の意味で良い推定量であると言えます。
これを平均対数尤度と呼んだりします。


もし、ここで新しいデータとして今の分布と同様の分布から Zが得られたとすると、

平均対数尤度の推定値は \frac{1}{n}\sum log(f(z; \hat{\theta}))となるわけですが、 \hat{\theta}はデータXの下での最尤推定量だったので、


 \frac{1}{n}\sum log(f(X; \hat{\theta)}) ≧ \frac{1}{n}\sum log(f(z; \hat{\theta)})

が成り立ちます。つまり、よりよい推定量を考えるためには平均対数尤度を最大化しなければならないのに、それよりも値が大きいXの対数尤度を最大化して満足してしまっていることがわかります。


しかし、そうはいっても「もしZを取ってこれたら」という話をしているに過ぎないので、実際はデータXでなんとかしなければなりません。

そこで、うまく平均対数尤度に近づくように補正する項を持ってきてやろうといことで漸近的に出てくるのが、先ほどAICで出てきた dim(\theta)です。

こうして、平均対数尤度を考える代わりに、既存のデータでどうにかできるようAICが生まれてきたというわけですね。

企業でのデータ解析は、「動くコードが書ける」だけではダメみたい

こんにちは!
最近ブログの更新が滞ってしまっていて申し訳ないです。

前にもお伝えしていた通り、某社で1か月ほどデータサイエンティストとしてインターンをさせていただいています。

今回は、そこで勤めてみて、初めて知ったことを雑記していこうと思います。


スポンサーリンク



界隈、機械学習屋さんが多い

まず、これについては就活をする中でも感じていましたが、データサイエンティスト界隈、機械学習屋さんが多いですね....

私は統計屋さんとして、インターンに参加したのですが、インターン生もほとんど機械学習屋、社内もほとんどが機械学習屋さんで、思っていたよりずっと比率が高い印象でした。

個人的には、言うて6:4くらいだろうと思っていたのですが、そんなもんじゃなさそうです。

まあ機械学習便利ですもんね。


引継ぎを意識して分析を整理しておく必要がある

これも、実際にお仕事をしてみて初めて知ったことなのですが、(というかまあ考えてみればそらそうだって感じですが)

”分析して、プレゼンして、面白い結果が出たので、実装しましょう終わり”

ではないんですね。その後、本人が居なくとも別の実装に活かせるように、引継ぎ資料を整理しておく必要があるようで、

今回のインターンで初めてmarkdownを書きました!一通りの分析を、書いてあるコードを実行することで試せるように、引継ぎ者に伝わるようにするためには資料作成の力もないといけないらしい。もう少しmarkdown資料、綺麗に作れるようになりたいなあ。


SQLも書ければ良いというわけでもない

SQL弄ってもらうから勉強してきてね~と事前に言われていたので、SQLもそれっぽく勉強していったのですが、

「うちの社内では、〇〇で規約を統一してるから、〇〇で動くコードをお願いね~」

とのこと。これも引継ぎのためでもあるのですが、やはりプロダクトにしやすくするために書き方も統一してるみたいで、その辺、自分勝手にコードを書いて、動けば良いじゃダメなんだ~!ってのも学びでした。プログラミング弱者の僕にとっては中々つらい事実ですが笑


分析関連の話や、データに関する話は多分言っちゃダメなのでその辺は控えておきますが、とにかく色々学びが多いです!
何かあまり書く時間もないので、雑な感じになっちゃいましたが、とりあえずこんな感じで!

また4月からはバリバリ書いてくんでよろしくお願いします~

そういえば、全然関係ないけど約束のネバーランド面白いですね!インターン終わりの本屋で手に取り、そこから一気に読んでしまいました!


最新刊

1巻

Rのplotを使ってNealの漏斗を可視化してみる

今日は、「StanとRでベイズ統計モデリング」という本の10章にあるNealの漏斗問題の可視化を簡単なコードだけで行ってみたいと思います。

ちなみにStanとRでベイズ統計モデリングとは、こちらの本です。

こちらの本、全てのコードがしっかりサポートページに掲載されていて、とても勉強になるのですが、10章にあるNealの漏斗の図は省略されていて(私が見た限り載っていなくて/2019年1月3日現在)、どうやって可視化すればいいんだろうと思ったのでやってみました。
とりあえず、綺麗な図というよりは確認出来ればいいので簡単なplotで試してみたいと思います。

スポンサーリンク


Nealの漏斗問題

まず、Nealの漏斗問題とは何かというと、ベイズモデルを作った際に、パラメータ同士対数事後周辺尤度が次のような漏斗状の形になることで生じるサンプリング問題を指します。
f:id:bananarian:20190103152542p:plain

このような状況になると、愚直にStanでモデリングを行っても、サンプルを効率よく取ってくることが出来ません。そこで再パラメータ化という、一旦独立な分布からサンプリングしてやったのち、変数変換を行う処置をとるわけです。

Rコード

では、コードを打っていきます。

まず、rstanの読み込み

library(rstan)

次にモデルを作っておきます

model="
parameters{
	real a;
	vector[1000] r;
}

model{
	a~normal(0,3);
	r~normal(0,exp(a/2));
}"

fit_a=stan(model_code=model,seed=1234,chains=1,iter=4000)
a1=rstan::extract(fit_a)$a
r1=rstan::extract(fit_a)$r[,1]

この状態でのサンプリングはあまりうまくいきません。それをplotでチェックします。

r=seq(-10,10,0.025)
a=seq(-10,10,0.025)
n=length(r)
lp1=matrix(,n,n)
for(i in r){
	for(j in a){
		lp1[which(r==i),which(a==j)]=log(dnorm(i,0,3))+log(dnorm(j,0,exp(i/2)))
	}
}
lp1[lp1==-Inf]=0
contour(r,a,t(lp1))
points(r1,a1,col="red")

f:id:bananarian:20190103153136p:plain

あまりうまくサンプリング出来ていないことが確認出来ました。

さて、再パラメータ化を行った場合はどうでしょう。

model2="
parameters{
	real a_raw;
	vector[1000] r_raw;
}

transformed parameters{
	real a;
	vector[1000] r;
	a=3.0*a_raw;
	r=exp(a/2)*r_raw;
}

model{
	a_raw~normal(0,1);
	r_raw~normal(0,1);
}"

fit_b=stan(model_code=model2,seed=1234,chains=1,iter=4000)
a2=rstan::extract(fit_b)$a
r2=rstan::extract(fit_b)$r[,1]

重ねてみる

contour(r,a,t(lp1))
points(r2,a2,col="green")

f:id:bananarian:20190103153354p:plain

漏斗の先も取ってきて来れていることが確認出来ました!

次元削減とは

今日は統計学や機械学習、その他画像認識等の場面で登場する「次元削減」の話について書こうと思います。


まず、次元削減とは、 n次元データを n-k次元にうつすことで、次元を下げることを指します。

簡単に例を出すと、
例えば国語・数学・理科・社会の4科目(1科目0点~100点)のテストをしたとします。
この時太郎君は国語:70点、数学:90点、理科:85点、社会:60点だったとします。

今現在変数は4つあって、
国語の得点(0点~100点)
数学の得点(0点~100点)
理科の得点(0点~100点)
社会の得点(0点~100点)となってるわけですが、

太郎君の点数を見た塾の先生や、太郎君の親は

「太郎君は、理系科目が得意だね~」

なんて言ったりしますよね。

スポンサーリンク



でも、今ある変数に「理系科目変数」なんてものは無いわけです。
当たり前ですが、これは次のようなプロセスで導かれています。

国語・社会⇒文系科目
数学・理科⇒理系科目

と置き換え、
(文系科目)<(理系科目)という何かしらの基準に基づく大小関係を考え、「太郎君は理系科目の方が得意」という話になるわけです。

特に数学や国語に重きを置くわけでないのなら、例えば次のような変換が出来ますね。
 X_{国語},X_{社会},X_{数学},X_{理科}という得点変数について、

 Y_{文系科目}=\frac{1}{2}X_{国語}+\frac{1}{2}X_{社会}
 Y_{理系科目}=\frac{1}{2}X_{数学}+\frac{1}{2}X_{理科}

今回であれば太郎君の Y_{文系科目},Y_{理系科目}は次のようになるはずです。
 Y_{文系科目} =65
 Y_{理系科目} =87.5

つまり、 Y_{文系科目}< Y_{理系科目}で、太郎君は理系科目の方が得意だとなるわけです。


この簡単な変換も実は「次元削減」です。

4次元(国語・数学・理科・社会)あったデータが2次元(文系科目・理系科目)になっていますよね。


この次元削減、どんな利点があるのでしょう。
何の利点も無いのに、人間が普段の会話で何となくこんな次元削減したりすることはないはずです。


一番の理由は多次元空間の解釈の難しさにあります。
国語・数学・理科・社会ならまださておき、8科目、9科目、かなり細かい専門に分けて10科目以上のテストをしたとして、その人の特性をどう解釈したらよいでしょうか。生データのままでは、せいぜい「この科目よりこの科目の方が点数が高い」とかそれくらいしか言えません。

結局色々なデータがあっても、多次元すぎてどう解釈すればわかりませんねーとなるわけです。

そこで、有用な情報を残したまま(あまり差の出ない情報は捨てて)、解釈しやすい次元までデータを持って行ってやろう

となるわけですね。これが次元削減です。


しかし、有用な情報を残すとは具体的にどういうことでしょう。
先ほどの例は、国語と社会は文系科目という分類が一般認識としてあるので、次元削減の方法までは考える必要もありませんでしたが、実データでは、「削減したいのはやまやまだけど、どう分けたらいいか分からない!」という状況の方が多かったりします。

そこで普通は(色々方法はありますが、よくある手法として)、データの分散に着目するわけですね。
主成分分析や決定木等の手法は、分散に着目した分析法です。


何故分散に着目するかというと、話は単純で

分散が小さい⇒データ間に決定的に差が無い

と考えられるからです。
例えば10人生徒に対して、国語のテストをしたとします。点数が次のようになったとすると
 10,11,23,21,5,90,89,95,99,100
国語の得意な生徒と苦手な生徒に分けるのは簡単ですね。つまりこのデータには、国語の得意・苦手の情報がしっかり入っているため、国語が得意か苦手か分けることが容易ということになります。

点数が次のようになっていたとすると
 76,70,69,71,77,83,65,72,72
これで国語の得意な生徒と苦手な生徒に分けるのはしんどそうです。
つまり、このデータは分散が小さく、皆が似通った点数を取っているため、国語の得意な生徒と苦手な生徒に分けるのは難しいということになります。

もしかしたら、生徒皆が国語のテストが得意で、全国平均20点のテストを受けているのかもしれませんし、実は生徒皆が国語のテストが苦手で、全国平均90点のテストを受けてこの点数になっているかもしれません。また、70点くらいまでならだれでも取れるけど、それ以上を取るのが難しいようなテストで、得意な人も苦手な人も点数が均されている可能性もありますね。

何が言いたいかというと、分散が小さい分類基準よりも、分散が大きい分類基準の方が、そのデータの性質が見出されやすいということです。

だから、特に主成分分析では、分散が大きくなるように主成分を決定していきますよね。

これが、次元削減の概要になります。
また次回にでも主成分分析等々の手法について記事をあげていこうと思います。

C++を使ってディープラーニングの基本を試す

最近、Rだけで統計解析したりすることに限界を感じているので、そうだ!C++を勉強しようと思い立ちました。

そこでとりあえず下の本を参考に、ディープラーニングの実装しながらC++を勉強しようなんていう次第です。

コードをそのまま写経するのも勉強にはなるのですが、私がやりたいのはディープラーニングを学ぶことではなく、
C++を学ぶことにあるので、少し修正する形でコードを組みなおしてみました。

スポンサーリンク



単純パーセプトロン

まず簡単に単純パーセプトロンの説明を。
 t_i=0 1かを分類する問題を考えます。

ある信号(ヒント?) X_iが与えられたら、出力される t_iを考えたいわけです。

で、何個か教師データを与えてやって学習させます。
例えば X_1=(1,0,0,1)の時は t_1=0になりましたよーみたいなデータを事前にたくさん学習させ、分類器を考えます。

何のモデルも無く分類器を考えることは出来ないので次のようなモデルを仮定します。

 W_i・X_i≧0ならば z_i=1
 W_i・X_i<0ならば z_i=0
 W_i=(w_0,w_1,...w_s)
ただし、 X_iの第一項は1


という分類器を作るわけですね。 zが予測値です。

昔、サポートベクターマシンのハードマージンに関する記事で線形分離可能の話をしましたが、この単純パーセプトロンでも線形分離可能を仮定します。
【初心者向け】2値分類SVM(サポートベクターマシン)の考え方と最適化問題 - バナナでもわかる話

で、実際の正解値と予測値の差が0になるような W_iを求めます。
線形分離可能なのでそういう組が見つかるはずですね。

C++コード

ということでコードを本を参考に少々修正しながら書いてみました。
ビジュアルスタジオのコンソールを使っています。

とりあえず名前空間stdを使うことにします。

#include "pch.h"
#include <iostream>
using namespace std;

parcepクラスを定義して単純パーセプトロンを行うクラスを設定してみました。
まだクラスを置くことの利点はよくわかってませんが、とりあえずC++の練習と言うことでw

class parcep {
public:
		float dot(float *v1, float *v2, int len);
		float step(float v);
		float forward(float *x, float *w, int len);
		void train(float *w, float *x, float t, float e, int len);
};

//2種類のベクトルの内積をとる関数
float parcep::dot(float *v1, float *v2, int len) {
	float sum = 0;
	for (int i = 0; i < len; i++) {
		sum += v1[i] * v2[i];
	}
	return sum;
}

//内積が0を超えたら1,超えなかったら0を出す関数
float parcep::step(float v) {
	return v > 0 ? 1 : 0;
}

//内積を取って、0か1かを出力する関数
float parcep::forward(float *x, float *w, int len) {
	float u = parcep::dot(x, w, len);
	return parcep::step(u);
}

//正解データtとデータによる出力zを比べて差がなくなるまでwをチューニングするトレーニング関数
void parcep::train(float *w, float *x, float t, float e, int len) {
	float z = parcep::forward(x, w, len);
	for (int j = 0; j < len; j++) {
		w[j] += (t - z)*x[j] * e;
	}
}

//データ数DATA_NUMSの定義
//ウエイトの個数WEIGH_NUMSの定義
#define DATA_NUMS 5
#define WEIGH_NUMS 4

//試してみる
	int main()

{
		parcep pa;
		//学習率を定義
		float e = 0.1;

		//入力データ
		float x[DATA_NUMS][WEIGH_NUMS]
			= { {1,0,0,1},{1,0,1,1},{1,0,1,0},{1,1,1,1},{1,0,0,0}};
		//正解データ
		float t[DATA_NUMS]
			= { 0,0,0,1,0 };
		//重みの初期値
		float w[WEIGH_NUMS]
			= { 0,0,0,0 };

		//トレーニング回数の設定
		int time = 10;
		for (int i = 0; i < time; i++) {
			cout << "time" << i;
			for (int j = 0; j < DATA_NUMS; j++) {
				pa.train(w, x[j], t[j], e, WEIGH_NUMS);
			}
			for (int j = 0; j < WEIGH_NUMS; j++) {
				cout <<endl<< "w" << j << ":" << w[j];
			}
			cout << endl;
		}

		for (int i = 0; i < DATA_NUMS; i++) {
			cout << pa.forward(x[i], w, WEIGH_NUMS);
		}
}


出力はこんな感じ
正解データが00010なので、訓練の末、パラメータ wが収束して、予測値も00010になったら正解です。

time0
w0:0
w1:0.1
w2:0.1
w3:0.1

time1
w0:-0.1
w1:0.1
w2:0.1
w3:0

time2
w0:-0.1
w1:0.1
w2:0.1
w3:0

time3
w0:-0.1
w1:0.1
w2:0.1
w3:0

time4
w0:-0.1
w1:0.1
w2:0.1
w3:0

time5
w0:-0.1
w1:0.1
w2:0.1
w3:0

time6
w0:-0.1
w1:0.1
w2:0.1
w3:0

time7
w0:-0.1
w1:0.1
w2:0.1
w3:0

time8
w0:-0.1
w1:0.1
w2:0.1
w3:0

time9
w0:-0.1
w1:0.1
w2:0.1
w3:0

00010

実際の予測データも00010となってうまくいっています。
パラメータも w0=-0.1,w1=0.1,w2=0.1,w3=0で止まっていますね。

画像認識始めました~Rでスクラッチでエッジ検出をやってみた~

ここ最近画像認識を勉強中です。
とりあえず読んだ傍からスクラッチで書いていっています!
参考本はこれ。

そして過去記事一覧はこちら
画像認識 カテゴリーの記事一覧 - バナナでもわかる話



前回はバイラテラルフィルタをやってみました。


今回はバイラテラルフィルタを使ってノイズを取った画像に対して、エッジ検出というものを行ってみます!

スポンサーリンク



エッジ検出

物体と物体の間の縁(エッジ)を検出する技術をエッジ検出と呼びます!
基本的には画像における画素差の急激な変化であったり輝度の急激な変化が起こった場所を物体の縁だと判断しています!
ただ、これもかなりノイズに引っ張られるので、前回前々回とやったフィルタリング技術を使ったり、検出時点でノイズを考慮したモデルを使う必要があるわけですね。

元画像

ちなみに無加工の元画像はこれです
f:id:bananarian:20181120002450p:plain

で、パイラテラルフィルタをかけたものがこれです。
f:id:bananarian:20181123081624p:plain

前者よりも後者の方がノイズが減っているのが分かりますね。後者を使っていきます。


微分

画素変化は微分で判断します。
ただ、画素値は離散個の点なので、よくある連続関数に対する微分は使えません。
ただ、要は1個ずれた時の差を検出すれば良いわけなので横方向、縦方向に関して差を取ってやれば良いということになります。

微分フィルタ

とりあえずx方向に関して微分するフィルタを用意してやります。

#x方向に対する微分フィルタFを作る
Fx=(1/2)*cbind(c(0,-1,0),rep(0,3),c(0,1,0))
F=array(Fx,dim=c(3,3,3))
#適用してみる
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image6[i,j,]=apply((F*zzz),3,sum)
	}
}

出力結果
f:id:bananarian:20181123085938p:plain

パソコンで見るとうっすら白い線でエッジが見えるのですが、アップしたらただ真っ黒な画像になってますね(笑)
画像のノイズがまだまだ多いので、ほぼエッジが検出出来ていないことによります。そこで、ノイズを考慮した手法を次に紹介します。


ちなみにy方向への微分はこんな感じ

Fy=(1/2)*cbind(rep(0,3),c(-1,0,1),rep(0,3))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image6[i,j,]=apply((F*zzz),3,sum)
	}
}

ノイズも考慮してエッジを検出する方法にプリューウィッドフィルタがあります。

プリューウィッドフィルタ

周りの画素の微分も考慮して周りのノイズの影響を減らそうという発想で作ったのがこれです。

#プリューウィッドフィルタ
image7=array(0,dim=c(dim(image)[1],dim(image)[2],3))
#x方向への微分
Fx=cbind(rep(-1,3),rep(0,3),rep(1,3))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image7[i,j,]=apply((F*zzz),3,sum)
	}
}

出力結果
f:id:bananarian:20181123100443p:plain


さっきよりも縁が浮き出てますね。手の形とかもわかるようになってます。
実際の画像だともっとキレイに縁検出出来てるんですけど、はてなブログに投稿すると見えにくくなりますね~。
投稿される時に圧縮されてるんでしょうか。

ちなみにy方向へのフィルタはこんな感じ

Fy=cbind(c(-1,0,1),c(-1,0,1),c(-1,0,1))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image7[i,j,]=apply((F*zzz),3,sum)
	}
}

スポンサーリンク



ソーベルフィルタ

プリューウィッドフィルタだと、周りと中心が均等に重みづけされて微分しているわけですが、当然中心が一番重要なのは間違いないので、中心に多めに重みを付けたフィルタがソーベルフィルタです。

コードを見てもらうとわかる通り、真ん中の部分が1じゃなくて2になってますね。

image8=array(0,dim=c(dim(image)[1],dim(image)[2],3))
#x方向への微分
Fx=cbind(c(-1,-2,-1),rep(0,3),c(1,2,1))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image8[i,j,]=apply((F*zzz),3,sum)
	}
}

出力結果
f:id:bananarian:20181123101410p:plain

プリューウィッドフィルタよりも縁が濃く検出されていますね。

ちなみにy方向はこんな感じ

Fy=cbind(c(-1,0,1),c(-2,0,2),c(-1,0,1))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image8[i,j,]=apply((F*zzz),3,sum)
	}
}

ラプラシアンフィルタ

ここまでは1階微分の話でしたが、二階微分して輝度で縁を検出することも出来ます。

image9=array(0,dim=c(dim(image)[1],dim(image)[2],3))
F0=cbind(c(0,1,0),c(1,-4,1),c(0,1,0))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image9[i,j,]=apply((F*zzz),3,sum)
	}
}

出力結果
f:id:bananarian:20181123102821p:plain

良い感じですね~

斜めも考慮したモデルを考えることも出来ます。

image9=array(0,dim=c(dim(image)[1],dim(image)[2],3))
F0=cbind(c(1,1,1),c(1,-8,1),c(1,1,1))
F=array(Fx,dim=c(3,3,3))
for(i in 2:(dim(image)[1]-1)){
	for(j in 2:(dim(image)[2]-1)){
		zzz=image5[max(1,(i-1)):min(dim(image)[1],(i+1)),max(1,(j-1)):min(dim(image)[2],(j+1)),]
		image9[i,j,]=apply((F*zzz),3,sum)
	}
}

f:id:bananarian:20181123103546p:plain


こんな感じで、画像の中にある物体のエッジを検出することが出来るわけですね。

画像認識始めました~バイラテラルフィルタをRで実装してみた

最近画像認識について勉強中です。とりあえず、数式で理解した後スクラッチで書いてみるのが一番早いかなと思ったので、読んだ傍からRで書いています。
画像認識始めました~手始めにRでスクラッチでフィルタリングしてみた~ - バナナでもわかる話
画像認識始めました~ガウシアンフィルタをRで実装してみた - バナナでもわかる話

参考本はこれ。

前回はガウシアンフィルタを使って前に撮ったイグアナちゃんをフィルタリングしたわけなんですけど、
今回はバイラテラルフィルタでやってみます。

スポンサーリンク


バイラテラルフィルタ

ガウシアンだと、各画素に対して均等に均していくことになるわけですが、これだと境界領域も均してしまうことになります。

そこで、画素の差にも考慮して重みづけしてやろうというのがこのフィルタです。


早速Rでコードを書いてみました。

#k個周りのバイラテラルフィルタ
k=40
#位置の分散
sigma=3
#画素の分散
sigmar=3
#imageに画像を挿入する(前記事参照)
#あとは実行するだけ
image5=array(0,dim=c(dim(image)[1],dim(image)[2],3))
bira=function(mat=matrix(),ii,jj,kk,sigma,sigmar){
	mid_gaso=mat[i,j,]
	zzz=mat[max(1,(ii-kk)):min(dim(image)[1],(ii+kk)),max(1,(jj-kk)):min(dim(image)[2],(jj+kk)),]
	hajix=max(1,(ii-kk))
	hajiy=max(1,(jj-kk))
	idox=1-hajix
	idoy=1-hajiy
	midPointx=ii+idox
	midPointy=jj+idoy
	Dim_matx=matrix(1:dim(zzz)[1],dim(zzz)[1],dim(zzz)[2])-midPointx
	Dim_maty=t(matrix(1:dim(zzz)[2],dim(zzz)[2],dim(zzz)[1]))-midPointy
	ggg=matrix(0,dim(zzz)[1],dim(zzz)[2])
	for(iii in 1:dim(zzz)[1]){
		for(jjj in 1:dim(zzz)[2]){
			ggg[iii,jjj]=exp(-(Dim_matx[iii,jjj]^2+Dim_maty[iii,jjj]^2)/(2*sigma^2)-sum((zzz[iii,jjj,]-mid_gaso)^2)/(2*sigmar^2))
		}
	}
	ggg
}

for(i in 1:(dim(image)[1])){
	for(j in 1:(dim(image)[2])){
		GG0=bira(mat=image,ii=i,jj=j,kk=k,sigma=sigma,sigmar=sigmar)
		GG=array(GG0,dim=c(dim(GG0)[1],dim(GG0)[2],3))
		zzz=image[max(1,(i-k)):min(dim(image)[1],(i+k)),max(1,(j-k)):min(dim(image)[2],(j+k)),]
		W=sum(GG0)
		image5[i,j,]=(1/W)*apply((GG*zzz),3,sum)
	}
}

出力結果
f:id:bananarian:20181123081624p:plain


もとのやつ
f:id:bananarian:20181120002450p:plain


かなりノイズが減って滑らかになっているのがわかるでしょうか。
こうやってノイズを消すんですね~。

次回はエッジの検出をやっていきます。

画像認識始めました~ガウシアンフィルタをRで実装してみた

最近画像認識について勉強中です。
画像認識始めました~手始めにRでスクラッチでフィルタリングしてみた~ - バナナでもわかる話


とりあえず、数式で理解した後スクラッチで書いてみるのが一番早いかなと思ったので、読んだ傍からRで書いています。
参考本はこれ。

前回は平均化フィルタを使って前に撮ったイグアナちゃんをフィルタリングしたわけなんですけど、
今回はガウシアンフィルタでやってみます。

スポンサーリンク



ガウシアンフィルタとは

二変量(画像には縦と横があるので)標準正規分布の中心の点(0,0)を平滑化したい点に合わせて、そこから離れるほど重みが小さくなるようなウエイトを用意してやろうという方法で、離散個のピクセルにウエイトを与えてウエイトの和が1にならねばならないので、その調整のために確率の総和 Wで割ります。条件付き分布と発想は同じですね。


というわけでコードを書いてみました。正直自分コーディング技術がないので、半ばゴリ押しのようになってしまいましたが....

#k個周りのガウシアンフィルタ
#パラメータの設定(sigmaは標準正規分布の標準偏差)
k=5
sigma=3
#imageに処理したい画像を入れてる
image4=array(0,dim=c(dim(image)[1],dim(image)[2],3))
#matに処理したい画像を入れる
Gaussian=function(mat=matrix(),ii,jj,kk,sigma){
	zzz=mat[max(1,(ii-kk)):min(dim(image)[1],(ii+kk)),max(1,(jj-kk)):min(dim(image)[2],(jj+kk)),]
#取り出した行列のスポットしたい点の座標を見つける
	hajix=max(1,(ii-kk))
	hajiy=max(1,(jj-kk))
	idox=1-hajix
	idoy=1-hajiy
#座標の特定
	midPointx=ii+idox
	midPointy=jj+idoy
	Dim_matx=matrix(1:dim(zzz)[1],dim(zzz)[1],dim(zzz)[2])-midPointx
	Dim_maty=t(matrix(1:dim(zzz)[2],dim(zzz)[2],dim(zzz)[1]))-midPointy
#gggにガウシアンフィルタを作る
	ggg=matrix(0,dim(zzz)[1],dim(zzz)[2])
	for(iii in 1:dim(zzz)[1]){
		for(jjj in 1:dim(zzz)[2]){
			ggg[iii,jjj]=(1/(2*pi*sigma^2))*exp(-(Dim_matx[iii,jjj]^2+Dim_maty[iii,jjj]^2)/(2*sigma^2))
		}
	}
	ggg
}
for(i in 1:(dim(image)[1])){
	for(j in 1:(dim(image)[2])){
		GG0=Gaussian(mat=image,ii=i,jj=j,kk=k,sigma=sigma)
		GG=array(GG0,dim=c(dim(GG0)[1],dim(GG0)[2],3))
		zzz=image[max(1,(i-k)):min(dim(image)[1],(i+k)),max(1,(j-k)):min(dim(image)[2],(j+k)),]
#Wはウエイトの和
		W=sum(GG0)
		image4[i,j,]=(1/W)*apply((GG*zzz),3,sum)
	}
}

こんな感じで作ってみました。コーディング技術欲しい.....

というわけで処理結果
f:id:bananarian:20181120122403p:plain

ちなみに処理前
f:id:bananarian:20181120002450p:plain

口の周りとか足のあたりとか見てもらうとわかりやすいかと思うんですけど、ザラザラ感がちょっと消えてます。
もう少しkを大きくしたりsigmaを大きく取ったりするともう少し分かりやすく変化するとは思いますが、


要はフィルタリングって前回も言ったと思うんですけど、ノイズを取ることが目的です。
だからザラザラした部分とかが均されてるわけです。

ただこのガウシアンフィルタの問題点は何か特徴に従って区別して均しているわけではなく、全ての点で同じような均し方をしているため、物体間の境界が不鮮明になる可能性があります。それを解決するためのフィルタがあるようなので、次回はそれを実装しようと思います。