作業：P64，2.5節，4：
考慮上機題2中的3個類別，設p（wi）=1/3 （a）以下各測試點與上機練習2中各類別均值間的Mahalanobis距離分別是多少： （1，2，1）t，（5，3，2）t，（0，0，0）t，（1，0，0）t （b）對以上各點進行分類。 （c）若設p（w1）=0.8，p（w2）=p（w3）=0.1，再對以上測試點進行分類。
這道題的思路完全與P35的例1一樣，所以會做例1，這道題就會做了。
我以以例題1的數據，測試了[0,0]，[3,8]分類正確。
另外，我不會python的= =，然后，我做完后才知道可以用java(╬▔皿▔)凸，這代碼是憋出來的，至今為止依然不會用自帶的矩陣乘法、轉置balabala，只能自己實現……。
這兩道題花了我二天半的時間（中途心態崩潰導致買買買導致錯過deadline，我正在寫博客的現在……qwq還在等老師的聯系方式……老天保佑）。
事實證明，如果只為了趕快交作業的話，就算你沒去上課，你也不用看書。 去做例1把，那里什么都有，如果遇到不懂的，就順着它的思路往前去找，比如說書才買的我，看不懂例1，就從第二章開始找u、sigma的定義，再盲猜例1中u、sigma的計算方式，才看懂例1。
（最新補充，我發現了，協方差矩陣的算法在P72第二行，注意x、u均為列向量）
（a）：
1、用式（40）計算出均值向量
2、用式（41）計算出協方差矩陣，
3、遍歷測試點和樣本數據
4、實現並獲得（x-u）、（x-u）.T
5、代入sqrt( (x-μ)'Σ^(-1)(x-μ) ) 計算Mahalanobis距離
（b）：
1、輸入樣本數據，得到u、sigma矩陣
2、遍歷測試點，對每個測試點，計算其分別屬於類i的判別函數值，選出最大的，判定為類i（判別函數方程對應《模式分類第二版》P32 （66）~（69））
（注意樣本數據和測試數據的格式和書上保持一直）
1 # a.py 2 import numpy 3 4 5 # 輸入：x0[i]第i個測試點，x[i]第i個樣本數據，group類總數，len0維度（描述一個類的元素總數） 6 # 輸出：dis[i][j]表示第i個測試點和第j個類的均值向量的距離 7 def mahalanobis(x0, x, group, len0): 8 dis = list() 9 u = getU(x, group, len0) 10 11 for i in range(len(x0)): 12 xx = numpy.array(x0[i]) 13 sigma = getSigma(x, group, len0) 14 dis.append([]) 15 for j in range(len(u)): 16 u[j] = numpy.array(u[j]) 17 18 # sqrt( (x-μ)'Σ^(-1)(x-μ) ) 計算Mahalanobis距離 19 t = [xx - u[j]] 20 t1 = myMult(t, numpy.linalg.inv(sigma[j])) 21 t2 = myMult(t1, rowToCol(xx - u[j])) 22 23 dis[i].append(t2[0][0]) 24 25 return dis 26 27 28 # 輸入：（同上） 29 # 輸出：u[i]是類i的均值向量 30 def getU(x, group, len0): 31 u = list() 32 for i in range(group): 33 t = list() 34 for j in range(len(x)): 35 t.append([]) 36 for z in range(i * len0, i * len0 + len0): 37 t[j].append(x[j][z]) 38 u.append(numpy.mean(t, 0)) 39 return u 40 41 42 # 輸入：（同上） 43 # 輸出：返回值sigma[i]表示第i類的協方差矩陣 44 def getSigma(x, group, len0): 45 sigma = list() 46 u = getU(x, group, len0) 47 48 for i in range(group): 49 sigma.append([]) 50 51 for i in range(group): 52 # t對所有樣本數據，只保留類i的列 53 t = list() 54 for j in range(len(x)): 55 t.append([]) 56 for z in range(i * len0, i * len0 + len0): 57 t[j].append(x[j][z]) 58 t[j] = t[j] - u[i] 59 sigma[i] = covariance(t) 60 return sigma 61 62 63 # 輸入：x對所有樣本數據，只保留某個類的列 64 # 輸出：某個類的協方差矩陣 65 def covariance(x): 66 all = covariance0(x[len(x) - 1]) 67 for i in range(len(x) - 1): 68 add(covariance0(x[i]), all) 69 70 for i in range(len(all)): 71 for j in range(len(all[0])): 72 all[i][j] /= len(x) 73 74 return all 75 76 # 輸入：a一個單獨的樣本數據，並且只保留某個類的列 77 # 輸出：a的協方差矩陣 78 def covariance0(a): 79 aa = [] 80 for i in range(len(a)): 81 aa.append([a[i]]) 82 83 ans = [] 84 for i in range(len(aa)): 85 t = [] 86 for j in range(len(aa)): 87 t.append(aa[i][0] * a[j]) 88 ans.append(t) 89 return ans 90 91 92 # ---------------------工具方法------------------------------ 93 94 95 # 行向量轉置為列向量 96 def rowToCol(x): 97 t = list(); 98 for i in range(len(x)): 99 t.append([]) 100 101 for i in range(len(x)): 102 t[i].append(x[i]) 103 return t 104 105 106 # 利用副作用，把矩陣a加到矩陣all里 107 def add(a, all): 108 for i in range(len(a)): 109 for j in range(len(a[0])): 110 all[i][j] += a[i][j] 111 112 113 # 矩陣乘法 114 def myMult(a, b): 115 ans = list() 116 for i in range(len(a)): 117 ans.append([]) 118 119 for i in range(len(a)): 120 for j in range(len(b[0])): 121 ii = 0 122 jj = 0 123 ans[i].append(0) 124 while ii < len(a[0]): 125 ans[i][j] = ans[i][j] + (a[i][ii] * b[jj][j]) 126 ii = ii + 1 127 jj = jj + 1 128 return ans 129 130 131 # 樣本數據 132 x = [ 133 numpy.array([-5.01, -8.12, -3.68, -0.91, -0.18, -0.05, 5.35, 2.26, 8.13]), 134 numpy.array([-5.43, -3.48, -3.54, 1.30, -2.06, -3.53, 5.12, 3.22, -2.66]), 135 numpy.array([1.08, -5.52, 1.66, -7.75, -4.54, -0.95, -1.34, -5.31, -9.87]), 136 numpy.array([0.86, -3.78, -4.11, -5.47, 0.50, 3.92, 4.48, 3.42, 5.19]), 137 numpy.array([-2.67, 0.63, 7.39, 6.14, 5.72, -4.85, 7.11, 2.39, 9.12]), 138 numpy.array([4.94, 3.29, 2.08, 3.60, 1.26, 4.36, 7.17, 4.33, -0.98]), 139 numpy.array([-2.51, 2.09, -2.59, 5.37, -4.63, -3.65, 5.75, 3.97, 6.65]), 140 numpy.array([-2.25, -2.13, -6.94, 7.18, 1.46, -6.66, 0.77, 0.27, 2.41]), 141 numpy.array([5.56, 2.86, -2.26, -7.39, 1.17, 6.30, 0.90, -0.43, -8.71]), 142 numpy.array([1.03, -3.33, 4.33, -7.50, -6.32, -0.31, 3.52, -0.36, 6.43]), 143 ] 144 145 # 測試點 146 x0 = [ 147 numpy.array([1, 2, 1]), 148 numpy.array([5, 3, 2]), 149 numpy.array([0, 0, 0]), 150 numpy.array([1, 0, 0]), 151 ] 152 153 ans = mahalanobis(x0, x, 3, 3) 154 155 for i in range(len(ans)): 156 print(ans[i])
1 # b.py 2 import numpy 3 4 5 # 輸入：x0[i]第i個測試點，x[i]樣本i，p[i]類i的先驗概率，group類別總數，len0維度 6 # 輸出：cla[i]第i個測試點的類別 7 def classify(x0, x, p, group, len0): 8 sigma = getSigma(x, group, len0) 9 u = getU(x, group, len0) 10 cla = list() 11 12 for i in range(len(x0)): 13 cla.append([]) 14 15 for i in range(0, len(x0)): 16 max = -100000000 17 # 根據最大的函數判斷值歸類 18 # 判別函數方程對應《模式分類第二版》P32 （66）~（69） 19 for g0 in range(0, group): 20 gg = g(x0[i], u[g0], sigma[g0], p[g0]) 21 if gg > max: 22 cla[i] = g0 + 1 23 max = gg 24 25 return cla 26 27 28 # 輸入：x測試點，u某個類的均值向量，sigma某個類的協方差矩陣，p某個類的先驗概率 29 # 輸出：x的判別函數值（方程對應《模式分類第二版》P32 （66）~（69）） 30 def g(x, u, sigma, p): 31 Wi = myMult2(-0.5, numpy.linalg.inv(sigma)) 32 t1 = myMult([x], Wi) 33 t2 = myMult(t1, rowToCol(x - u))[0][0] 34 35 wi = myMult(numpy.linalg.inv(sigma), rowToCol(u)) 36 t3 = myMult([colToRow(wi)], rowToCol(x))[0][0] 37 38 t4 = myMult([u], numpy.linalg.inv(sigma)) 39 t5 = myMult(t4, rowToCol(u)) 40 t6 = myMult2(-0.5, t5)[0][0] 41 t6 = t6 - 0.5 * numpy.log(numpy.linalg.det(sigma)) + numpy.log(p) 42 43 return t2 + t3 + t6 44 45 46 # ---------------------下列函數已經出現在a.py中---------------------------- 47 48 def getU(x, group, len0): 49 u = list() 50 for i in range(group): 51 t = list() 52 for j in range(len(x)): 53 t.append([]) 54 for z in range(i * len0, i * len0 + len0): 55 # print(j, z) 56 t[j].append(x[j][z]) 57 u.append(numpy.mean(t, 0)) 58 return u 59 60 61 def getSigma(x, group, len0): 62 sigma = list() 63 u = getU(x, group, len0) 64 65 for i in range(group): 66 sigma.append([]) 67 68 for i in range(group): 69 t = list() 70 for j in range(len(x)): 71 t.append([]) 72 for z in range(i * len0, i * len0 + len0): 73 t[j].append(x[j][z]) 74 t[j] = t[j] - u[i] 75 sigma[i] = covariance(t) 76 return sigma 77 78 79 def covariance(x): 80 all = covariance0(x[len(x) - 1]) 81 for i in range(len(x) - 1): 82 add(covariance0(x[i]), all) 83 84 for i in range(len(all)): 85 for j in range(len(all[0])): 86 all[i][j] /= len(x) 87 88 return all 89 90 91 def covariance0(a): 92 aa = [] 93 for i in range(len(a)): 94 aa.append([a[i]]) 95 96 aa 97 ans = [] 98 for i in range(len(aa)): 99 t = [] 100 for j in range(len(aa)): 101 t.append(aa[i][0] * a[j]) 102 ans.append(t) 103 return ans 104 105 106 # ------------------工具方法------------------ 107 108 109 # 利用副作用，把a加到all里 110 def add(a, all): 111 for i in range(len(a)): 112 for j in range(len(a[0])): 113 all[i][j] += a[i][j] 114 115 116 # 行向量轉置為列向量 117 def rowToCol(x): 118 t = list(); 119 for i in range(len(x)): 120 t.append([]) 121 122 for i in range(len(x)): 123 t[i].append(x[i]) 124 return t 125 126 127 # 列向量轉置為行向量 128 def colToRow(x): 129 t = list() 130 for i in range(len(x)): 131 t.append(x[i][0]) 132 return t 133 134 135 # 矩陣乘法 136 def myMult(a, b): 137 ans = list() 138 for i in range(len(a)): 139 ans.append([]) 140 141 for i in range(len(a)): 142 for j in range(len(b[0])): 143 ii = 0 144 jj = 0 145 ans[i].append(0) 146 while ii < len(a[0]): 147 ans[i][j] = ans[i][j] + (a[i][ii] * b[jj][j]) 148 ii = ii + 1 149 jj = jj + 1 150 return ans 151 152 153 # 常數*矩陣 154 def myMult2(n, x): 155 ans = list() 156 for i in range(len(x)): 157 ans.append([]) 158 for j in range(len(x[0])): 159 ans[i].append([]) 160 ans[i][j] = x[i][j] * n 161 return ans 162 163 164 a = [ 165 numpy.array([-5.01, -8.12, -3.68, -0.91, -0.18, -0.05, 5.35, 2.26, 8.13]), 166 numpy.array([-5.43, -3.48, -3.54, 1.30, -2.06, -3.53, 5.12, 3.22, -2.66]), 167 numpy.array([1.08, -5.52, 1.66, -7.75, -4.54, -0.95, -1.34, -5.31, -9.87]), 168 numpy.array([0.86, -3.78, -4.11, -5.47, 0.50, 3.92, 4.48, 3.42, 5.19]), 169 numpy.array([-2.67, 0.63, 7.39, 6.14, 5.72, -4.85, 7.11, 2.39, 9.12]), 170 numpy.array([4.94, 3.29, 2.08, 3.60, 1.26, 4.36, 7.17, 4.33, -0.98]), 171 numpy.array([-2.51, 2.09, -2.59, 5.37, -4.63, -3.65, 5.75, 3.97, 6.65]), 172 numpy.array([-2.25, -2.13, -6.94, 7.18, 1.46, -6.66, 0.77, 0.27, 2.41]), 173 numpy.array([5.56, 2.86, -2.26, -7.39, 1.17, 6.30, 0.90, -0.43, -8.71]), 174 numpy.array([1.03, -3.33, 4.33, -7.50, -6.32, -0.31, 3.52, -0.36, 6.43]), 175 ] 176 177 a0 = [ 178 numpy.array([1, 2, 1]), 179 numpy.array([5, 3, 2]), 180 numpy.array([0, 0, 0]), 181 numpy.array([1, 0, 0]), 182 ] 183 184 # p = numpy.array([1 / 3, 1 / 3, 1 / 3]) 185 p = numpy.array([0.8, 0.1, 0.1]) 186 187 c = classify(a0, a, p, 3, 3) 188 print(c)