本实验包括三个部分,一是求解最大似然估计的均值和方差,二是Parzen 窗估计和设计分类器及KNN概率密度估计,三是KNN方法来构建一个分类器。以上三个实验的数据及完整代码:https://github.com/JIANG54864/ml_exp/tree/main/exp2

一、 最大似然估计

【实验内容及说明】

类1 类2
样本 x1 x2 x3 x1 x2 x3
1 0.76 -2.1 -1.9 0.49 0.0012 -0.065
2 -1 -2.6 0.87 0.054 -0.063 -0.12
3 -5.3 -0.15 -1.6 0.47 0.034 -0.35
4 0.89 -4.7 -0.029 0.69 0.1 0.17
5 1.9 2.2 -0.23 0.55 -0.18 -0.011
6 -0.3 -0.87 0.27 0.61 0.12 -0.27
7 -0.087 0.58 0.42 0.58 0.089 -0.4
8 -3.3 -3.4 -0.2 0.27 -0.04 -0.31
9 -0.32 1.7 1.3 0.055 -0.035 0.38
10 0.71 0.23 0.39 0.53 0.011 -0.15

使用上面给出的三维数据:
(1)编写程序,对类 1 和类 2 中的三个特征分别求解最大似然估计的均值和方差。
(2)编写程序,处理二维数据的情形。对类 1 和类 2 中任意两个特征的组合分别求解最大似然估计的均值和方差(每个类有3种可能)。
(3)编写程序,处理三维数据的情形$p(x)\sim N(\mu,\Sigma)$,对类 1 和类 2 中三个特征求解最大似然估计的均值和方差。
(4)假设该三维高斯模型是可分离的,编写程序估计类 1 和类 2 中的均值和协方差矩阵中的参数。
(5)比较前 4 种方法计算出来的每一个特征的均值的异同,并加以解释。
(6)比较前 4 种方法计算出来的每一个特征的方差的异同,并加以解释。

【正文】

概率论中学过,最大似然估计是一种参数估计的方法,利用已知样本结果,反推最大可能导致这样结果的参数值,有了这些参数值就可以用来做贝叶斯决策或者其它任务。最大似然估计的均值和方差定义和代码计算如下(多维情况,一维情况简单计算均值和方差即可):

$ \bar\mu = \frac{1}{n}\sum_{k=1}^{n}\bar{x_k} $

故均值可用以下函数计算:

1
2
3
4
5
avg = pd.Series(index=vectors.columns, dtype=float)
for column in vectors.columns:
# 计算每个特征的平均值
avg[column] = vectors[column].mean()
return np.array(avg)

$ \bar\Sigma=\frac{1}{n}\sum_{k=1}^{n}{(\bar{x_k}-\bar\mu)(\bar{x_k}-\bar\mu)^T} $

方差可用以下函数计算:

1
2
3
4
5
6
7
8
9
def calculateCov(vectors: pd.DataFrame):
mu = np.matrix(calculateAvg(vectors)).T
dimension = vectors.shape[1]
Cov = np.zeros((dimension, dimension))
for index, row in vectors.iterrows():
xi = np.matrix(row).T
diff = xi - mu
Cov += diff * diff.T
return Cov / vectors.shape[0]

三维高斯模型是可分离的,即当三维高斯模型的协方差矩阵是对角矩阵时,各特征间相互独立,联合概率分布可分解为各特征边缘分布的乘积。此时,协方差矩阵的非对角线元素(协方差)为零,对角线元素为各特征的方差。这意味着各维度之间没有相关性,模型参数估计时仅需考虑各维度自身的均值和方差。

关于均值和方差的比较,均值相同,所有方法均采用样本均值,与协方差结构无关。方差相同,差异仅体现在协方差项(非对角线元素)的存在与否,方差估计为样本方差,仅依赖各维度数据,与协方差矩阵是否包含非零协方差无关。

二、非参数估计

【实验内容及说明】

类1 类2 类3
样本 x1 x2 x3 x1 x2 x3 x1 x2 x3
1 -1.63 1.31 0.28 1.03 -2.18 0.08 0.68 0.14 2.17
2 -0.78 0.58 1.52 1.96 1.27 0.45 1.41 0.92 0.13
3 5.73 0.21 -0.81 3.12 0.13 0.16 1.22 0.69 0.99
4 -4.32 1.18 -0.44 1.23 -0.21 -0.11 2.46 -0.38 1.45
5 -6.2 0.82 0.07 1.39 -1.44 -0.19 0.85 0.87 0.79
6 2.77 3.16 0.65 0.82 0.34 -0.16 2.51 1.35 3.22
7 4.38 2.42 2.20 0.94 -1.38 -0.21 0.60 1.31 2.44
8 6.21 2.01 0.91 1.28 -0.12 0.17 0.64 0.97 2.19
9 -0.96 1.93 1.54 2.31 0.011 0.14 1.36 0.99 0.58
10 -0.19 1.94 -0.26 1.94 0.26 0.08 0.66 0.88 0.51

Parzen 窗估计:
使用上面表格中的数据进行 Parzen 窗估计和设计分类器。窗函数为一个球形的高斯函数如下所示:
img
编写程序,使用 Parzen 窗估计方法对任意一个的测试样本点img进行分类。对分类器的训练则使用表格中的三维数据。令h = 1,分类样本点为$(0.5,1.0,0.0)^𝑡,(0.31,1.51,-0.50)^𝑡,(-0.3,0.44, -0.1)^𝑡$ 。

k-近邻概率密度估计:
对上面表格中的数据使用k-近邻方法进行概率密度估计:
1) 编写程序,对于一维的情况,当有 n 个数据样本点时,进行k-近邻概率密度估计。对表格中的类3的特征𝑥1,用程序画出当 k=1,3,5 时的概率密度估计结果。
2) 编写程序,对于二维的情况,当有 n 个数据样本点时,进行k-近邻概率密度估计。对表格中的类2的特征$(𝑥1, 𝑥2)^𝑡$,用程序画出当 k=1,3,5 时的概率密度估计结果。
3) 编写程序,对表格中的3个类别的三维特征,使用k-近邻概率密度估计方法。并且对下列点处的概率密度进行估计:$ (-0.41,0.82,0.88)^𝑡,(0.14,0.72, 4.1)^𝑡,(-0.81,0.61, -0.38)^𝑡$。

【正文】

在参数估计中,我们总是假设概率密度函数的参数形式已知,但在现实世界中我们一般给出的概率密度的形式很少符合实际情况。

“非参数化方法”(non- parametric method)能处理任意的概率分布而不必假设密度的参数形式已知。

在模式识别中有多种令人感兴趣的非参数化方法。其中之一是如何从训练样本中估计概率密度函数$p(x|w_j)$。如果这种估计的结果是可靠的话那么在设计分类器时估计出的结果就可以认为是真正的概率密度。另一种方法讨论如何直接估计后验概率$P(w_j|x)$。这种方法的实现方式与所设计的算法直接相关,比如说“最近邻规则”就省略了概率估计这一步而直接进行判别函数的设计。

对于Parzen 窗估计,窗函数已给出,直接计算即可:

1
2
3
4
5
def window(sample: pd.Series, trainSample: pd.Series, h):
vector_s = np.matrix(sample).T# 转置
vector_ts = np.matrix(trainSample).T
diff = vector_s - vector_ts # 差值
return math.exp(-diff.T * diff / (2 * h ** 2)) # 按公式算

遍历训练集中的每一行,调用window函数计算测试样本与训练样本之间的窗函数值,并累加到似然值中,计算平均似然值并返回:

1
2
3
4
5
6
def Parzen(sample: pd.Series, trainSet: pd.DataFrame):
likelihood = 0.0
for index, row in trainSet.iterrows():
likelihood += window(sample, row, 1)
likelihood = likelihood / len(trainSet) #取平均值
return likelihood

计算样本在三个不同训练集上的条件概率密度。比较三个条件概率密度,输出样本属于概率密度最大的类别即可完成分类。

使用k-近邻方法进行概率密度估计时可以选择欧式距离或者马氏距离,欧氏距离适合特征独立、数据分布均匀、计算资源有限的场景。马氏距离适合特征相关性强、数据分布复杂、对准确性要求较高的场景。

欧式距离:$d(x,y)=\sqrt{\sum\limits_{i=1}^n(x_i-y_i)^2}$

马氏距离:$d(x,y)=\sqrt{(x-y)^TS^{-1}(x-y)}$,其中S是协方差矩阵

观察上面给出的数据,可以得出分布是比较均匀的,并且为了简化计算,就先使用欧氏距离了。

一维:

1
2
3
4
5
6
7
8
9
10
def one_dimension_KNN(testData: float, trainSet: pd.Series, k: int):
distance = []
for i in range(trainSet.shape[0]): # 遍历训练集,shape[0]表示数据的行数
d = np.abs(testData - trainSet[i])
distance.append(d)
distance.sort()
posterior = (k / trainSet.shape[0]) / (2 * distance[k - 1])
# 分子 k / trainSet.shape[0] 表示标点周围最近的 k 个样本,占训练集总样本的比例。分母 2 * distance[k - 1] 表示目标点附近的区间长度。
# 相除得到概率密度
return posterior

二维:

1
2
3
4
5
6
7
8
9
10
11
12
13
def two_dimension_KNN(x1: np.matrix, x2: np.matrix, trainSet: pd.DataFrame, k: int):
posteriorMatrix = np.zeros((x1.shape[0], x1.shape[1]))
for i in range(x1.shape[0]):
for j in range(x1.shape[1]):
distance = []
for index, row in trainSet.iterrows():
d = np.sqrt((x1[i, j] - row[0]) ** 2 + (x2[i, j] - row[1]) ** 2)
distance.append(d)
distance.sort()
posterior = (k / trainSet.shape[0]) / (np.pi * (distance[k - 1] ** 2) + np.spacing(1))
# 同样是计算距离,只不过是二维欧氏距离
posteriorMatrix[i, j] = posterior
return posteriorMatrix

三维:

1
2
3
4
5
6
7
8
9
10
11
12
13
def three_dimension_KNN(testData: np.matrix, k: int):
distance = [[], [], []]
posterior = []
for i in range(len(trainSet)):
for j in range(trainSet[i].shape[0]):
# 仍然是计算欧氏距离,只不过是三维
d = np.sqrt((testData[0, 0] - trainSet[i].iloc[j]['x1']) ** 2 +
(testData[1, 0] - trainSet[i].iloc[j]['x2']) ** 2 +
(testData[2, 0] - trainSet[i].iloc[j]['x3']) ** 2)
distance[i].append(d)
distance[i].sort()
V = 4 * np.pi * (distance[i][k - 1] ** 3) / 3 # 计算球体体积
posterior.append(k / (trainSet[i].shape[0]) / V)

分类时,对每个类求条件概率密度,最大的即为数据点预测的类

三、KNN实战

现有一诊断数据集diagnosis_result.csv,共有100条数据,数据格式如下图所示:

img

id代表序列,diagnosis_result表示Negative/Positive(阴性或阳性),需要通过后面8个属性特征来判断诊断结果为阴性或阳性。用学过的KNN方法来构建一个分类器,判断一个样本的诊断结果,KNN 模型需要自己实现,不可使用现成的第三方库

【正文】

首先进行归一化处理,为了后续接受新数据进行预测,还要保存各特征的最大最小值以便对新数据进行归一化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 归一化,先获取最大值最小值
# 获取除了 'id' 和 'result' 之外的特征名称
feature_keys = [col for col in reader.fieldnames if col not in ['id', 'diagnosis_result']]
# print(feature_keys)
max_min_values = {} # 存储各特征的最大最小值

for key in feature_keys:
values = [float(data[key]) for data in datas]
max_val = max(values)
min_val = min(values)
max_min_values[key] = (max_val, min_val) # 保存到字典
# 归一化处理
for data in datas:
data[key] = (float(data[key]) - min_val) / (max_val - min_val)

计算测试数据与训练集中每个数据的距离,按距离从小到大排序,并取前K个最近的数据。
对这K个数据进行加权平均,距离越小权重越大。根据加权结果返回分类结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def knn(data):
res = [
{"result": train["diagnosis_result"], "distance": distance(data, train)}
for train in train_set
]
res = sorted(res, key=lambda item: item['distance'])
res2 = res[0:k]
result = {'Positive': 0, 'Negative': 0}
sum_dist = 0
for r1 in res2:
sum_dist += r1['distance']
# 逐个分类加和
for r2 in res2:
result[r2["result"]] += 1 - (r2["distance"] / sum_dist)
if result['Positive'] > result['Negative']:
return 'Positive'
else:
return 'Negative'

对新数据进行预测:

1
2
3
4
5
6
7
def predict_new_data(new_data):
normalized_data = {}
for key in feature_keys:
max_val, min_val = max_min_values[key]
normalized_data[key] = (float(new_data[key]) - min_val) / (max_val - min_val)
normalized_data['diagnosis_result'] = 'Unknown'
return knn(normalized_data)

参考文献:

[1] Richard O. Duda, Peter E. Hart, David G. Stork 著;李宏东,等 译. 模式分类(原书第2版)[M]. 北京: 机械工业出版社, 2003.