作业:P64,2.5节,4:

 

考虑上机题2中的3个类别,设p(wi)=1/3

(a)以下各测试点与上机练习2中各类别均值间的Mahalanobis距离分别是多少:
(121)t,(532)t,(000)t,(100)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)