ロジスティック回帰（多クラス分類）について

#ロジスティック回帰識別器をpythonで自力で作りたい

２クラスのロジスティック回帰識別器をnumpyを用いて作成する事ができました。
次に、３クラス以上に分類する為にシグモイド関数ではなくソフトマックス関数を使い、それに関する誤差の微分則などの数学的な仕組みは理解済みです。

pythonでの実装で困っていること。

入力をx=(x_1,x_2,...,x_D)^Tとします。（D×1ベクトル）とし、出力がK個（クラスの数だけ）になるようなソフトマックス関数の実装方法がわからず困っています。（数式は理解しています。pythonでの実装方法がいまいちぴんと来ていません）

今回実装したいソフトマックスの式の形

（追記）
画像中のw、bは識別に使う重み（パラメータ）です。

2クラスでの識別器（自分で実装したもの）

python
1import numpy as np
2import makeGaussianData
3import matplotlib.pyplot as plt
4
5K = 2
6w = 0.02*np.random.rand(K)-0.01 #1.パラメータの初期化
7b=0.02*np.random.rand()-0.01
8
9X, lab, t = makeGaussianData.getData(K)
10z = np.empty(X.shape[0],)
11h = np.empty(X.shape[0],)
12eta=0.01
13cnt = 0
14
15def s(x):
16    return 1/(1+np.exp(-x))
17
18for i in range(10000): #2.適当な回数の繰り返し（学習データは400個）
19     n = np.random.randint(0,400) #i.400個のデータからランダムで1つ選択
20     z[n]= s(w @ X[n] + b)
21     h[n]= (-t[n]*np.log(z[n]))-((1-t[n])*np.log(1-z[n])) #ii.モデルの出力を求める
22     w = w-(eta*(X[n]*(z[n]-t[n])))
23     b = b-(eta*(z[n]-t[n])) #iii.パラメータの更新
24     if i % 1000 == 0: #2.iv. 1000の倍数になったときの処理
25        for j in range(X.shape[0]):
26            z[j]=s(w @ X[j] + b)
27            h[j]= (-1*t[j])*np.log(z[j])-(1-t[j])*np.log(1-z[j])
28            if z[j] > 0.5:
29                T = 1
30            else:
31                T = 0
32            if T == t[j]:
33                cnt = cnt+1
34
35        H = np.mean(h)
36        print("#{0}, H:{1} , {2}/{3}={4}".format(i,H,cnt,X.shape[0],cnt/X.shape[0]))
37        cnt=0
38
39fig = plt.figure()
40plt.xlim(-0.2, 1.2)
41plt.ylim(-0.2, 1.2)
42ax = fig.add_subplot(1, 1, 1)
43ax.set_aspect(1)
44ax.scatter(X[lab == 0, 0], X[lab == 0, 1], color = 'red')
45ax.scatter(X[lab == 1, 0], X[lab == 1, 1], color = 'green')
46if K == 3:
47    ax.scatter(X[lab == 2, 0], X[lab == 2, 1], color = 'blue')
48
49fig.show()
50plt.show()

makeGaussianData.py

python
1import numpy as np
2
3
4def getData(nclass, seed = None):
5
6    assert nclass == 2 or nclass == 3
7
8    if seed != None:
9        np.random.seed(seed)
10
11    # 2次元の spherical な正規分布3つからデータを生成
12    X0   = 0.10 * np.random.randn(200, 2) + [ 0.3, 0.3 ]
13    X1   = 0.10 * np.random.randn(200, 2) + [ 0.7, 0.6 ]
14    X2   = 0.05 * np.random.randn(200, 2) + [ 0.3, 0.7 ]
15
16    # それらのラベル用のarray
17    lab0 = np.zeros(X0.shape[0], dtype = int)
18    lab1 = np.zeros(X1.shape[0], dtype = int) + 1
19    lab2 = np.zeros(X2.shape[0], dtype = int) + 2
20
21    # X （入力データ）, label （クラスラベル）, t（教師信号） をつくる
22    if nclass == 2:
23        X = np.vstack((X0, X1))
24        label = np.hstack((lab0, lab1))
25        t = np.zeros(X.shape[0])
26        t[label == 1] = 1.0
27    else:
28        X = np.vstack((X0, X1, X2))
29        label = np.hstack((lab0, lab1, lab2))
30        t = np.zeros((X.shape[0], nclass))
31        for ik in range(nclass):
32            t[label == ik, ik] = 1.0
33
34    return X, label, t
35    
36
37if __name__ == '__main__':
38
39    import matplotlib
40    import matplotlib.pyplot as plt
41
42    K = 3
43    
44    X, lab, t = getData(K)
45
46    fig = plt.figure()
47    plt.xlim(-0.2, 1.2)
48    plt.ylim(-0.2, 1.2)
49    ax = fig.add_subplot(1, 1, 1)
50    ax.set_aspect(1)
51    ax.scatter(X[lab == 0, 0], X[lab == 0, 1], color = 'red')
52    ax.scatter(X[lab == 1, 0], X[lab == 1, 1], color = 'green')
53    if K == 3:
54        ax.scatter(X[lab == 2, 0], X[lab == 2, 1], color = 'blue')
55    plt.show()