作业：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)