サポートベクターマシンのパーセプトロンアルゴリズムをPythonで
2クラス分類を解く一番簡単なアルゴリズムです。
#適当な10点のx1,x2成分 x1 = np.array([0,2,4,6,8,10,12,14,16,18]) x2 = np.array([1,2,1,4,3,9,7,6,7,11]) #それぞれのラベル y = np.array([-1,-1,-1,-1,-1,1,1,1,1,1]) #点を全てプロット plt.scatter(px, py) #重みは0に近い値で初期化 w = np.array([0.001,0.001,0.001]) #失敗の数を表す k = 0 #サンプルの中での、ノルムの最大値 R = np.linalg.norm(np.array([x1[9], x2[9]])) #学習率 適当に小さく eta = 0.001 #繰り返し継続判定のフラグ miss = True #パーセプトロンの学習 while(miss): miss = False for i in range(10): if(y[i] * (w[0]*x1[i] + w[1]*x2[i] + w[2]) <= 0): #分類ミス miss = True w[0] = w[0] + eta * y[i] * x1[i] w[1] = w[1] + eta * y[i] * x2[i] w[2] = w[2] + eta * y[i] * R * R k += 1 #-3から20まで0.1刻みでnumpy配列で用意 x = np.arange(-3, 20, 0.1) #学習して得た直線の式 y = -(w[0] / w[1])*x - (w[2] / w[1]) plt.plot(x, y) plt.show()
いい感じに分類してくれていますね!
最小二乗法をPythonで実装
最小二乗法では、重みはこうなります。
\begin{align}
w\bf = (X^TX)^{-1}X^Ty\bf
\end{align}
#適当な6点のx,y座標 px = np.array([0,1,2,3,4,5]) py = np.array([1,5,6,8,9,9]) e = np.array([1,1,1,1,1,1]) #成分が全て1のベクトル X = np.array([px, e]).T #6点のx成分と1を縦に並べた行列(最小二乗法に必要) #最小二乗法による重みの解 w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(py) #点を全てプロット plt.scatter(px, py) #-3から7まで0.1刻みでnumpy配列で用意 x = np.arange(-3, 7, 0.1) #得られた重みを用いて、直線の式を作る y = w[0]*x + w[1] #得られた直線を引く plt.plot(x, y) plt.show()
ちなみに、Windrow-Hoffアルゴリズムという最小二乗法を解く手法も紹介する
#適当な6点のx,y座標 px = np.array([0,1,2,3,4,5]) py = np.array([1,5,6,8,9,9]) #重みは0で初期化 w = np.array([0, 0]) #学習率 適当に小さく eta = 0.001 #Widrow-Hoffアルゴリズムで最小二乗法を解く for _ in range(1000): for i in range(6): w = w - eta*(w[0]*px[i] + w[1] - py[i])*np.array([px[i], 1]) #点を全てプロット plt.scatter(px, py) #-3から7まで0.1刻みでnumpy配列で用意 x = np.arange(-3, 7, 0.1) #得られた重みを用いて、直線の式を作る y = w[0]*x + w[1] #得られた直線を引く plt.plot(x, y) plt.show()
点に沿ったいい感じの直線が引けてて気持ちいいね!
線形計画問題
線形計画問題の解法の一つである「シンプレックス法」をPythonで実装してみました。
def LP(c, A, b, xb_idx): xn_idx = np.arange(A.shape[1]) xn_idx = np.delete(xn_idx, xb_idx) while True: B = A[:, xb_idx] N = A[:, xn_idx] B_inv = np.linalg.inv(B) # Bの逆行列 b_bar = np.dot(B_inv, b) # bバー。(xB, xN) = (bバー, 0) if np.any(b_bar < 0) : print('initial basic solution isn\'t solution!') return 0 cB = c[xb_idx] cN = c[xn_idx] pi = np.dot(B_inv.T, cB) # シンプレックス乗数 rcf = cN - np.dot(N.T, pi) # 相対コスト係数 if np.all(rcf >= 0) : solution = np.zeros_like(c, dtype=float) solution[xb_idx] = b_bar return solution for i in range(len(rcf)): if rcf[i] < 0 : choose = xn_idx[i] choose_idx = i break y = np.dot(B_inv, A[:, choose]) y = y.reshape(len(y), 1) if np.all(y <= 0) : print('Unbounded') return 0 theta_opt = b_bar / y theta = np.min(theta_opt) dropout = np.argmin(theta_opt) temp = xn_idx[choose_idx] xn_idx[choose_idx] = xb_idx[dropout] xb_idx[dropout] = temp xb_idx = sorted(xb_idx) xn_idx = sorted(xn_idx)
具体的にはこう使います。
\begin{align}
&minimize &-2x_1-3x_2 \\
&\text{subject to} &4x_1+3x_2+x_3&=120\\
&&x_1+2x_2+x_4&=60\\
&&x_i \geq 0 &&\forall{i = 1, 2, 3, 4}
\end{align}
こんな時は
A = np.array([[4,3,1,0],[1,2,0,1]]) b = np.array([[120,60]]).T c = np.array([[-2,-3,0,0]]).T xb_idx = [2,3] x = LP(c, A, b, xb_idx) print('Optimum solution') print(x) print('objective function -> ',np.round(np.dot(c.T, x)[0][0])) print('constraints -> ', np.all(np.round(np.dot(A, x) - b, 7) == 0))
出力は
Optimum solution
[[12.]
[24.]
[ 0.]
[ 0.]]
objective function -> -96.0
constraints -> True
matplotlibで陰関数のグラフを描画してみた
陰関数とは、というものです。
Pythonのmatplotlibを使って、を、陰関数に直して図示してみました。
まずは の定義から。
の元 に対して、そのを
\begin{align}
\|x\|_p = \sqrt[p]{|x_1|^p + |x_2|^p + \cdots + |x_n|^p}
\end{align}と定義する。
さて、定義から
\begin{align}
\|x\|_1 = |x_1| + |x_2| + \cdots + |x_n|
\end{align}
\begin{align}
\|x\|_2 = \sqrt{|x_1|² + |x_2|² + \cdots + |x_n|²}
\end{align}
ですね。
とし、の場合で図示してみましょう。
この場合、
import matplotlib.pyplot as plt import numpy as np #Lpノルムの p p = 200 x1range = np.arange(-2, 2, 0.025) x2range = np.arange(-2, 2, 0.025) X1, X2 = np.meshgrid(x1range, x2range) plt.axis([-2, 2, -2, 2]) plt.gca().set_aspect('equal', adjustable='box') Z = abs(X1)**p + abs(X2)**p - 1 plt.contour(X1, X2, Z, [0]) plt.show()
で、正方形っぽくなりましたね。
実は、で完全な正方形になります。
\begin{align}
\|x\|_\infty = max \left\{|x_1|, |x_2|, \cdots \right\}
\end{align}
の時、
なら、
なら、
となるため、
\begin{align}
\|x\|_\infty = 1 \longrightarrow |x_i| \leq 1 (i = 1, 2, \cdots, n)
\end{align}
となり、このようなグラフになるわけですね。
ではまた。