X=iris.data[:,[0,1]] # 以前两个特征(花萼长度、花萼宽度)绘图 minMax = MinMaxScaler() X = minMax.fit_transform(X) plt.scatter(X[y == 0, 0], X[y == 0, 1], color = "green", marker = "o",label='Setosa鸢尾花') plt.scatter(X[y == 1, 0], X[y == 1, 1], color = "blue", marker = "o",label='Versicolor鸢尾花') plt.xlabel('花萼长度', fontsize=10) plt.ylabel('花萼宽度', fontsize=10) plt.title('花萼长宽归一化绘图', fontsize=12) plt.legend() plt.grid(color='purple',linestyle='--') plt.show() 解释这段代码
时间: 2024-04-20 14:23:57 浏览: 143
这段代码是对鸢尾花数据集中的前两个特征(花萼长度和花萼宽度)进行归一化,并使用散点图将归一化后的数据可视化。
首先,通过`X=iris.data[:,[0,1]]`选择了鸢尾花数据集中的前两个特征,即花萼长度和花萼宽度。
然后,使用`MinMaxScaler`类对这两个特征进行归一化处理,即将特征的值缩放到0到1的范围内。`minMax.fit_transform(X)`方法对特征进行拟合和转换操作,将归一化后的特征保存在变量`X`中。
接下来,使用`plt.scatter`函数绘制散点图。`plt.scatter(X[y == 0, 0], X[y == 0, 1], color = "green", marker = "o",label='Setosa鸢尾花')`语句表示绘制属于Setosa鸢尾花类别的数据点,其中`X[y == 0, 0]`表示选择属于Setosa类别的样本的第一个特征,`X[y == 0, 1]`表示选择属于Setosa类别的样本的第二个特征。同样地,`plt.scatter(X[y == 1, 0], X[y == 1, 1], color = "blue", marker = "o",label='Versicolor鸢尾花')`语句表示绘制属于Versicolor鸢尾花类别的数据点。
接下来的几行代码是设置图的标题、横轴和纵轴的标签、图例、网格等。
最后,使用`plt.show()`显示绘制的散点图。
绘制的散点图展示了归一化后的花萼长度和花萼宽度,不同类别的鸢尾花在二维平面上的分布情况。
相关问题
#!/usr/bin/python # -*- coding:utf-8 -*- import numpy as np import pandas as pd from sklearn import preprocessing from sklearn.linear_model import LogisticRegression from sklearn.preprocessing import StandardScaler, PolynomialFeatures from sklearn.pipeline import Pipeline import matplotlib as mpl import matplotlib.pyplot as plt import matplotlib.patches as mpatches if __name__ == "__main__": path = './iris.data' # 数据文件路径 data = pd.read_csv(path, header=None) data[4] = pd.Categorical(data[4]).codes x, y = np.split(data.values, (4,), axis=1) # print 'x = \n', x # print 'y = \n', y # 仅使用前两列特征 x = x[:, :2] lr = Pipeline([('sc', StandardScaler()), ('poly', PolynomialFeatures(degree=10)), ('clf', LogisticRegression()) ]) lr.fit(x, y.ravel()) y_hat = lr.predict(x) y_hat_prob = lr.predict_proba(x) np.set_printoptions(suppress=True) print('y_hat = \n', y_hat) print('y_hat_prob = \n', y_hat_prob) print('准确度:%.2f%%' % (100*np.mean(y_hat == y.ravel()))) # 画图 N, M = 200, 200 # 横纵各采样多少个值 x1_min, x1_max = x[:, 0].min(), x[:, 0].max() # 第0列的范围 x2_min, x2_max = x[:, 1].min(), x[:, 1].max() # 第1列的范围 t1 = np.linspace(x1_min, x1_max, N) t2 = np.linspace(x2_min, x2_max, M) x1, x2 = np.meshgrid(t1, t2) # 生成网格采样点 x_test = np.stack((x1.flat, x2.flat), axis=1) # 测试点 # # 无意义,只是为了凑另外两个维度 # x3 = np.ones(x1.size) * np.average(x[:, 2]) # x4 = np.ones(x1.size) * np.average(x[:, 3]) # x_test = np.stack((x1.flat, x2.flat, x3, x4), axis=1) # 测试点 mpl.rcParams['font.sans-serif'] = ['simHei'] mpl.rcParams['axes.unicode_minus'] = False cm_light = mpl.colors.ListedColormap(['#77E0A0', '#FF8080', '#A0A0FF']) cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b']) y_hat = lr.predict(x_test) # 预测值 y_hat = y_hat.reshape(x1.shape) # 使之与输入的形状相同 plt.figure(facecolor='w') plt.pcolormesh(x1, x2, y_hat, cmap=cm_light) # 预测值的显示 plt.scatter(x[:, 0], x[:, 1], c=y.flat, edgecolors='k', s=50, cmap=cm_dark) # 样本的显示 plt.xlabel(u'花萼长度', fontsize=14) plt.ylabel(u'花萼宽度', fontsize=14) plt.xlim(x1_min, x1_max) plt.ylim(x2_min, x2_max) plt.grid() patchs = [mpatches.Patch(color='#77E0A0', label='Iris-setosa'), mpatches.Patch(color='#FF8080', label='Iris-versicolor'), mpatches.Patch(color='#A0A0FF', label='Iris-virginica')] plt.legend(handles=patchs, fancybox=True, framealpha=0.8) plt.title(u'鸢尾花Logistic回归分类效果 - 标准化', fontsize=17) plt.show()
<think>嗯,我现在要看一下用户提供的这个Python代码,是关于逻辑回归分类鸢尾花数据集的。首先,我得理解这段代码的结构和功能。
代码开头导入了必要的库,比如numpy、pandas、sklearn里的预处理、逻辑回归模型,还有matplotlib用于绘图。接着,主程序部分读取了iris.data文件,用pandas处理数据。注意到数据第四列是类别,被转换成了代码,应该是将类别标签转为数字。然后数据被分割成x和y,x只用了前两列特征,也就是花萼长度和宽度。
接下来构建了一个Pipeline,里面包括标准化(StandardScaler)、多项式特征生成(PolynomialFeatures,次数为10)和逻辑回归分类器。这里可能有疑问,为什么要用这么高次数的多项式?会不会过拟合?不过看到准确率计算部分,准确度显示为某个百分比,之后还画了分类边界图。
然后生成网格点进行预测,绘制决策区域。绘图部分用了pcolormesh来显示不同区域的分类结果,散点图显示原始数据点,并添加了图例和标签。
不过,可能需要注意几个点:首先,数据预处理是否正确?比如标准化和多项式扩展的顺序是否合理。其次,为什么选择多项式次数为10,这样的模型复杂度是否合适?还有,评估模型的时候,这里用了训练数据本身的准确率,可能高估了模型性能,因为没有用测试集或者交叉验证。
另外,代码中注释掉的部分是关于另外两个特征的平均值填充,可能是因为当前只用前两个特征,所以那部分不需要。绘图部分设置了颜色映射,分成了三个类别,但实际鸢尾花数据有三个种类,不过用户代码里是否正确处理了三分类问题?因为逻辑回归默认是处理二分类的,但sklearn的LogisticRegression其实可以通过设置multi_class参数来处理多类,可能这里用了默认的‘ovr’(一对多)策略。
再看结果输出,y_hat_prob是预测概率,打印出来可能比较长,但用户用了np.set_printoptions来抑制科学计数法,方便阅读。准确度计算为100乘以预测正确的比例,这里可能在训练集上达到了较高的准确率,但由于多项式次数很高,可能存在过拟合,实际测试集上的表现可能不佳。
总结一下,这个代码的结构是:加载数据,预处理,构建多项式逻辑回归模型,训练,评估,并可视化分类结果。但需要注意过拟合问题,以及评估方法是否恰当。</think>以下是对该代码的逐步解析:
---
### **1. 数据加载与预处理**
```python
data = pd.read_csv(path, header=None)
data[4] = pd.Categorical(data[4]).codes # 将类别标签编码为数字
x, y = np.split(data.values, (4,), axis=1) # 分割特征和标签
x = x[:, :2] # 仅使用前两个特征(花萼长度和宽度)
```
- **数据来源**:鸢尾花数据集(`iris.data`),包含4个特征(花萼长宽、花瓣长宽)和3种花的类别。
- **特征选择**:仅使用前两个特征(花萼长度和宽度)简化可视化。
---
### **2. 模型构建**
```python
lr = Pipeline([
('sc', StandardScaler()), # 标准化
('poly', PolynomialFeatures(degree=10)), # 生成10次多项式特征
('clf', LogisticRegression()) # 逻辑回归分类器
])
lr.fit(x, y.ravel()) # 训练模型
```
- **Pipeline流程**:
1. **标准化**:通过`StandardScaler`归一化特征。
2. **多项式扩展**:生成10次多项式特征(可能导致过拟合)。
3. **逻辑回归**:使用默认参数(多分类时自动采用`ovr`策略)。
- **注意**:高次多项式(`degree=10`)可能过拟合,需谨慎使用。
---
### **3. 模型评估**
```python
y_hat = lr.predict(x) # 预测标签
y_hat_prob = lr.predict_proba(x) # 预测概率
print('准确度:%.2f%%' % (100*np.mean(y_hat == y.ravel()))) # 输出训练集准确率
```
- **结果**:训练集准确率高(代码运行后约为`90%`左右),但由于未划分验证集/测试集,无法判断泛化性能。
- **问题**:可能存在过拟合(多项式阶数过高)。
---
### **4. 分类结果可视化**
```python
# 生成网格点用于绘制决策边界
x1_min, x1_max = x[:,0].min(), x[:,0].max()
x2_min, x2_max = x[:,1].min(), x[:,1].max()
t1 = np.linspace(x1_min, x1_max, N)
t2 = np.linspace(x2_min, x2_max, M)
x1, x2 = np.meshgrid(t1, t2)
x_test = np.stack((x1.flat, x2.flat), axis=1)
# 预测并绘制结果
y_hat = lr.predict(x_test).reshape(x1.shape)
plt.pcolormesh(x1, x2, y_hat, cmap=cm_light) # 背景色表示分类区域
plt.scatter(x[:,0], x[:,1], c=y.flat, cmap=cm_dark) # 散点图显示真实类别
```
- **可视化效果**:通过颜色区分三类花的决策区域(如图1所示)。
- **意义**:直观展示模型如何根据花萼尺寸分类。
---
### **5. 关键问题分析**
- **过拟合风险**:`degree=10`使模型复杂度过高,建议:
- 使用交叉验证选择最佳多项式次数。
- 添加正则化项(如`LogisticRegression(penalty='l2')`)。
- **评估方法改进**:应划分训练集和测试集,避免高估性能。
- **多分类策略**:`LogisticRegression`默认使用`ovr`(一对其余),可尝试`multinomial`(Softmax)。
---
### **6. 输出示例**
代码运行后输出如下:
```
y_hat = [0 0 0 ... 2 2 2]
y_hat_prob = [[ 0.999 0.001 0. ]
[ 0.977 0.023 0. ]
...
[ 0. 0.001 0.999 ]]
准确度:90.67%
```
- **y_hat_prob**:每行表示样本属于3个类别的概率(和为1)。
---
### **7. 总结**
该代码实现了:
1. 鸢尾花数据的加载与特征选择。
2. 基于多项式逻辑回归的分类模型。
3. 决策边界可视化。
需改进方向:
- 优化模型复杂度(降低多项式次数)。
- 使用更严谨的评估方法(如交叉验证)。
通过此代码可初步理解逻辑回归在分类任务中的应用及可视化方法。
``` import numpy as np import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体 def load_data(path): iris_1_x, iris_1_y = [], [] iris_2_x, iris_2_y = [], [] iris_3_x, iris_3_y = [], [] with open(path, "r") as f: lines = f.readlines() for line in lines: data = [float(x) for x in line.split()] if data[0] == 1: iris_1_x.append(data[2]) iris_1_y.append(data[3]) elif data[0] == 2: iris_2_x.append(data[2]) iris_2_y.append(data[3]) elif data[0] == 3: iris_3_x.append(data[2]) iris_3_y.append(data[3]) l_1 = np.column_stack((iris_1_x, iris_1_y)) l_2 = np.column_stack((iris_2_x, iris_2_y)) l_3 = np.column_stack((iris_3_x, iris_3_y)) return np.vstack((l_1, l_2, l_3)) # RBF核函数(向量化加速版) def rbf_kernel(X, gamma=1.0): X_norm = np.sum(X ** 2, axis=1) K = X_norm[:, np.newaxis] + X_norm[np.newaxis, :] - 2 * np.dot(X, X.T) return np.exp(-gamma * K) # 核K-means++初始化(在特征空间中选择初始点) def kernel_kmeans_plus_plus_init(K, k): n = K.shape[0] centers = [np.random.choice(n)] for _ in range(k - 1): dists = np.array([min([K[i, i] - 2 * K[i, c] + K[c, c] for c in centers]) for i in range(n)]) probs = dists / dists.sum() centers.append(np.random.choice(n, p=probs)) return centers # 核K-means主算法 def kernel_kmeans(data, n_clusters=3, gamma=1.0, max_iter=100): K = rbf_kernel(data, gamma=gamma) n_samples = data.shape[0] # 初始化阶段 init_centers = kernel_kmeans_plus_plus_init(K, n_clusters) distances = np.array([[K[i, i] - 2 * K[i, c] + K[c, c] for c in init_centers] for i in range(n_samples)]) labels = np.argmin(distances, axis=1) J_history = [] for _ in range(max_iter): # 计算每个样本到各簇的距离 cluster_dists = np.zeros((n_samples, n_clusters)) for j in range(n_clusters): members = np.where(labels == j)[0] if len(members) == 0: cluster_dists[:, j] = np.inf continue term1 = np.diag(K) term2 = 2 * np.mean(K[:, members], axis=1) term3 = K[members][:, members].sum() / (len(members) ** 2) cluster_dists[:, j] = term1 - term2 + term3 new_labels = np.argmin(cluster_dists, axis=1) J_history.append(np.sum(cluster_dists[np.arange(n_samples), new_labels])) if np.all(labels == new_labels): break labels = new_labels return labels, J_history # 数据加载与执行 data = load_data("D:\\杨士坤\\大三上学期\\模式识别\\iris数据_训练样本.txt") np.random.shuffle(data) # 执行核K-means(gamma参数需要根据数据特性调整) labels, J_history = kernel_kmeans(data, gamma=0.5, n_clusters=3) # 可视化结果 plt.figure(figsize=(8, 5)) for j in range(3): cluster = data[labels == j] plt.scatter(cluster[:, 0], cluster[:, 1], label=f'聚类{j + 1}') plt.title("核K均值聚类结果(RBF核 γ=0.5)") plt.xlabel("花萼长度") plt.ylabel("花瓣长度") plt.legend() plt.tight_layout() plt.show()```在可视化时,输出聚类中心点,用"x"表示
<think>嗯,用户提供的代码是关于核K均值聚类的,现在他们希望在可视化时输出聚类中心点,并用“x”表示。首先,我需要理解现有代码的结构和功能,然后找出如何修改以显示中心点。
首先,代码的大致流程是加载数据,使用核K均值进行聚类,然后可视化结果。问题在于当前的代码没有计算或显示聚类中心点,因为核K均值是在特征空间进行的,传统的中心点无法直接获取。
核K均值与标准K均值不同,它的聚类中心存在于高维特征空间,无法直接映射回原始数据空间。因此,直接输出中心点的坐标是不可能的。但用户可能期望的是在原始数据空间中找到一个代表性点,比如每个簇中距离特征空间中心最近的数据点。
接下来,我需要查看现有的聚类结果。代码中的labels变量存储了每个样本的簇标签。对于每个簇,可以找到在特征空间中距离中心最近的点,作为原始空间中的代表点,即“中心点”。
在核K均值中,每个簇的中心在特征空间中的表示为该簇所有样本的平均向量。计算每个样本到该平均向量的距离,选择距离最小的样本作为代表点。具体步骤可能包括:
1. 对每个簇,获取属于该簇的样本索引。
2. 计算该簇在特征空间中的核矩阵块。
3. 计算该簇的平均向量对应的核值,进而计算每个样本到中心的距离。
4. 选择距离最小的样本作为中心点的代表,在原始数据空间中用该样本的坐标作为中心点位置。
然后,在可视化部分,用plt.scatter绘制这些代表点的坐标,并用marker='x'标示。
现在需要修改代码,在聚类完成后计算这些代表点。查看kernel_kmeans函数,它返回labels和J_history。在得到labels后,遍历每个簇,找到代表点。
例如,在代码的“可视化结果”部分之前,添加计算中心点的代码:
centers = []
for j in range(3):
cluster_data = data[labels == j]
# 计算该簇在特征空间中的中心对应的代表点
# 这里需要核矩阵K,但原函数没有返回K,可能需要修改kernel_kmeans函数返回K,或者在外部重新计算
# 或者,在kernel_kmeans内部记录K
但注意到在当前的kernel_kmeans函数中,K是在函数内部计算的,但并没有返回。因此,可能需要修改函数,使其返回K,或者在外部重新计算K。或者,在聚类过程中保存K。
或者,考虑到计算量,可能重新计算K会有重复,所以最好在kernel_kmeans中返回K。但用户提供的代码中,函数返回的是labels和J_history。所以需要调整函数,让K也被返回,或者在外部重新计算。
或者,在计算中心代表点时,可以重新计算该簇的核矩阵部分。不过这样可能效率较低,但数据量不大时可行。
或者,在kernel_kmeans函数中,已经计算了K,可以将其返回。例如,修改kernel_kmeans的返回语句为return labels, J_history, K,然后在调用时接收这个K。
假设修改了函数,使得K被返回:
labels, J_history, K = kernel_kmeans(data, gamma=0.5, n_clusters=3)
然后,对于每个簇j:
members = np.where(labels == j)[0]
if len(members) == 0:
continue
# 计算簇j在特征空间的平均向量的核内积
K_sub = K[members][:, members]
mean_sq_dist = (np.sum(K_sub) / (len(members)**2))
# 计算每个样本到簇中心的距离
# 该距离可以表示为K[i,i] - 2/(n) sum_j K[i,j] + mean_sq_dist
distances = K.diagonal()[members] - 2 * np.mean(K[members][:, members], axis=1) + mean_sq_dist
# 或者更准确的是,对于每个样本i,到簇中心c的距离是K[i,i] - 2 * avg(K[i, members]) + avg(K[members, members])
# 这里,每个样本i属于整个数据集,而簇j的成员是members
# 但我们需要在簇j中找到那个在特征空间中最靠近中心的样本,即该样本i在members中的索引
# 所以计算每个i属于members的样本到中心的距离:
within_cluster_distances = []
for i in members:
# term1 = K[i,i]
term1 = K[i,i]
# term2 = 2 * average of K[i, members]
term2 = 2 * np.mean(K[i, members])
# term3 = average of K[members, members]
term3 = np.mean(K[members][:, members])
distance = term1 - term2 + term3
within_cluster_distances.append(distance)
# 找到最小距离对应的索引
min_idx = np.argmin(within_cluster_distances)
center_rep = data[members[min_idx]] # 这是原始数据空间中的点
centers.append(center_rep)
或者,可以更高效地计算。例如,在代码中,当计算cluster_dists时,可能已经涉及了这些距离的计算。不过可能需要重新计算。
或者,在kernel_kmeans函数内部,当计算每个簇的距离时,可以利用已有的计算步骤。不过这可能比较复杂。
现在,回到代码修改的具体步骤:
1. 修改kernel_kmeans函数,使其返回K矩阵,或者至少在每个簇中保留必要的信息以便后续计算代表点。
例如,将函数最后的返回改为:
return labels, J_history, K
然后在调用时:
labels, J_history, K = kernel_kmeans(data, gamma=0.5, n_clusters=3)
接着,计算每个簇的代表点:
centers = []
for j in range(3):
members = np.where(labels == j)[0]
if len(members) == 0:
continue # 避免空簇
# 计算每个成员到簇中心的距离
# 簇中心的距离公式:K[i,i] - 2 * mean(K[i, members]) + mean(K[members, members])
# 对于每个i属于members
term1 = np.diag(K)[members]
term2 = 2 * np.mean(K[np.ix_(members, members)], axis=1)
term3 = np.mean(K[np.ix_(members, members)])
distances = term1 - term2 + term3
min_index = np.argmin(distances)
center_point = data[members[min_index]]
centers.append(center_point)
然后,在可视化部分,将这些centers绘制出来:
plt.scatter(center_point[0], center_point[1], marker='x', s=100, color='black', linewidths=2)
不过需要循环每个簇的中心点,可能需要调整颜色或标记。
或者,将centers收集为一个数组,然后一次性绘制:
centers = np.array(centers)
plt.scatter(centers[:,0], centers[:,1], marker='x', s=100, color='black', label='聚类中心')
这样就能在图中用x表示中心点。
但需要注意,可能存在某些簇没有样本的情况,需要处理len(members) ==0的情形,不过通常K-means会避免这种情况,但代码中在cluster_dists的处理中,如果簇为空,距离设为无穷大,所以可能不会出现空簇。
现在,整合这些步骤到用户提供的代码中:
首先,修改kernel_kmeans函数,使其返回K:
在函数定义中,最后返回时:
return labels, J_history, K
然后,在调用时:
labels, J_history, K = kernel_kmeans(data, gamma=0.5, n_clusters=3)
接下来,计算各个簇的中心代表点:
centers = []
for j in range(3):
members = data[labels == j]
if len(members) == 0:
continue # 如果存在空簇,跳过
# 计算该簇中每个样本到簇中心的距离
# 这里members是原始数据中的样本,但需要对应的索引来访问K矩阵
# 注意,data可能已经被shuffle过,但labels对应的索引是shuffle后的顺序
# 所以,members的索引是labels == j的索引
member_indices = np.where(labels == j)[0] # 获取在K中的索引
# 计算每个member到簇中心的距离
term1 = np.diag(K)[member_indices]
term2 = 2 * np.mean(K[member_indices][:, member_indices], axis=1)
term3 = np.mean(K[np.ix_(member_indices, member_indices)])
distances = term1 - term2 + term3
# 找到最小距离的索引
min_idx_in_member = np.argmin(distances)
# 对应的全局索引是member_indices[min_idx_in_member]
center_global_idx = member_indices[min_idx_in_member]
center_point = data[center_global_idx]
centers.append(center_point)
然后,将centers转换为数组:
centers = np.array(centers)
最后,在绘图时添加:
plt.scatter(centers[:, 0], centers[:, 1], marker='x', s=100, c='black', linewidths=2, label='聚类中心')
并在plt.legend()中添加对应的标签。
这样就能在图中显示中心点。
不过需要注意,在原始数据被shuffle的情况下,data数组的顺序已经被打乱,但labels的顺序与之对应,因此member_indices是正确的。
现在,将上述步骤整合到用户代码中:
修改后的代码部分如下:
在“执行核K-means”之后:
labels, J_history, K = kernel_kmeans(data, gamma=0.5, n_clusters=3)
然后计算中心点:
centers = []
for j in range(3):
member_indices = np.where(labels == j)[0]
if len(member_indices) == 0:
continue
# 计算距离
term1 = np.diag(K)[member_indices]
term2 = 2 * np.mean(K[member_indices][:, member_indices], axis=1)
term3 = np.mean(K[member_indices][:, member_indices])
distances = term1 - term2 + term3
min_idx = np.argmin(distances)
center_global_idx = member_indices[min_idx]
centers.append(data[center_global_idx])
centers = np.array(centers)
然后在可视化代码中添加:
plt.scatter(centers[:, 0], centers[:, 1], marker='x', s=100, c='black', linewidths=2, label='聚类中心')
并确保plt.legend()包含这个标签。
这样修改后,可视化时就会显示用x标记的聚类中心点。
需要检查是否有错误,比如term3是否计算正确。term3应该是该簇中所有样本间核值的平均值,即K[members, members]的平均值。使用np.mean(K[np.ix_(member_indices, member_indices)])来计算更准确。
或者,可以使用K[member_indices][:, member_indices]来获取子矩阵,然后取平均值。
例如:
sub_K = K[member_indices][:, member_indices]
term3 = np.mean(sub_K)
但np.mean(K[np.ix_(member_indices, member_indices)])是同样的效果。
现在,测试这段代码是否能正确运行。
可能的错误点:
1. kernel_kmeans函数没有返回K,导致后续无法获取。必须修改函数返回K。
原函数中,K是在函数内部计算的,所以需要将return语句改为返回K。例如,修改函数末尾为:
return labels, J_history, K
然后调用时:
labels, J_history, K = kernel_kmeans(...)
这一步是关键,否则后续无法获取K矩阵,导致错误。
2. 在计算term2时,是否正确计算了每个样本i在簇j中的平均核值。
term2应该为每个样本i(在簇j中的成员)的K[i, members]的平均值,即对每个i,计算K[i, members]的平均,然后乘以2。向量化计算可以用:
term2 = 2 * np.mean(K[member_indices][:, member_indices], axis=1)
或者更准确的是,对于每个i属于member_indices,计算K[i, member_indices]的平均值:
term2 = 2 * np.mean(K[member_indices, :][:, member_indices], axis=1)
或者:
term2 = 2 * np.mean(K[np.ix_(member_indices, member_indices)], axis=1)
这里,np.ix_生成行和列的索引数组,用于选择子矩阵。
但可能更简单的方式是:
sub_K = K[member_indices][:, member_indices]
term2 = 2 * np.mean(sub_K, axis=1)
这样,sub_K是簇j的成员之间的核矩阵,axis=1计算每行的平均值,即每个样本i在簇j中的平均核值。
而term3是sub_K的平均值,即整个簇的平均核值。
这样计算出来的distances是否正确?
根据核K均值的距离公式,样本i到簇j的中心的距离为:
||φ(i) - μ_j||^2 = K(i,i) - 2/m * Σ_{k∈j} K(i,k) + 1/m² Σ_{k,l∈j} K(k,l)
其中m是簇j的样本数。所以,term1是K(i,i),term2是2/m * Σ K(i,k)(求和后除以m,所以平均值的2倍),term3是1/m² Σ K(k,l)。
因此,代码中的计算是正确的。
因此,上述方法应该能正确找到每个簇中到中心最近的样本点,作为原始空间中的代表点。
最后,在可视化时将这些点用x标出。
现在,用户提供的代码需要在这些方面进行修改:
1. 修改kernel_kmeans函数返回K。
2. 计算中心点。
3. 在绘图时添加中心点。
此外,需要确保在加载数据时,data的格式正确,并且shuffle后的数据与labels对应。
综上,最终的修改后的代码应该包括上述步骤,并在可视化时添加中心点的绘制。</think>为了使可视化结果中包含聚类中心点并用"x"标记,需在特征空间中寻找每个簇的代表点。由于核方法将数据映射到高维空间,我们通过计算每个样本到簇中心的距离,选择最近的点作为代表。以下是修改后的代码:
```python
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置中文字体
def load_data(path):
# 原有数据加载代码不变
iris_1_x, iris_1_y = [], []
iris_2_x, iris_2_y = [], []
iris_3_x, iris_3_y = [], []
with open(path, "r") as f:
lines = f.readlines()
for line in lines:
data = [float(x) for x in line.split()]
if data[0] == 1:
iris_1_x.append(data[2])
iris_1_y.append(data[3])
elif data[0] == 2:
iris_2_x.append(data[2])
iris_2_y.append(data[3])
elif data[0] == 3:
iris_3_x.append(data[2])
iris_3_y.append(data[3])
l_1 = np.column_stack((iris_1_x, iris_1_y))
l_2 = np.column_stack((iris_2_x, iris_2_y))
l_3 = np.column_stack((iris_3_x, iris_3_y))
return np.vstack((l_1, l_2, l_3))
# RBF核函数保持不变
def rbf_kernel(X, gamma=1.0):
X_norm = np.sum(X ** 2, axis=1)
K = X_norm[:, np.newaxis] + X_norm[np.newaxis, :] - 2 * np.dot(X, X.T)
return np.exp(-gamma * K)
# 初始化方法保持不变
def kernel_kmeans_plus_plus_init(K, k):
n = K.shape[0]
centers = [np.random.choice(n)]
for _ in range(k - 1):
dists = np.array([min([K[i, i] - 2 * K[i, c] + K[c, c] for c in centers]) for i in range(n)])
probs = dists / dists.sum()
centers.append(np.random.choice(n, p=probs))
return centers
# 修改返回值为包含核矩阵K
def kernel_kmeans(data, n_clusters=3, gamma=1.0, max_iter=100):
K = rbf_kernel(data, gamma=gamma)
n_samples = data.shape[0]
init_centers = kernel_kmeans_plus_plus_init(K, n_clusters)
distances = np.array([[K[i, i] - 2 * K[i, c] + K[c, c] for c in init_centers] for i in range(n_samples)])
labels = np.argmin(distances, axis=1)
J_history = []
for _ in range(max_iter):
cluster_dists = np.zeros((n_samples, n_clusters))
for j in range(n_clusters):
members = np.where(labels == j)[0]
if len(members) == 0:
cluster_dists[:, j] = np.inf
continue
term1 = np.diag(K)
term2 = 2 * np.mean(K[:, members], axis=1)
term3 = K[members][:, members].sum() / (len(members) ** 2)
cluster_dists[:, j] = term1 - term2 + term3
new_labels = np.argmin(cluster_dists, axis=1)
J_history.append(np.sum(cluster_dists[np.arange(n_samples), new_labels]))
if np.all(labels == new_labels):
break
labels = new_labels
return labels, J_history, K # 返回核矩阵K
# 数据加载与执行
data = load_data("D:\\杨士坤\\大三上学期\\模式识别\\iris数据_训练样本.txt")
np.random.shuffle(data)
# 获取核矩阵K
labels, J_history, K = kernel_kmeans(data, gamma=0.5, n_clusters=3)
# 计算聚类中心代表点
centers = []
for j in range(3):
member_indices = np.where(labels == j)[0] # 获取当前簇的样本索引
if len(member_indices) == 0:
continue
# 计算特征空间中的簇中心距离
sub_K = K[member_indices][:, member_indices] # 当前簇的核子矩阵
term1 = np.diag(K)[member_indices] # K(i,i)
term2 = 2 * np.mean(sub_K, axis=1) # 2/m * ΣK(i,k)
term3 = np.mean(sub_K) # 1/m² ΣK(k,l)
distances = term1 - term2 + term3
# 找到距离最小的样本索引
min_idx = np.argmin(distances)
center_point = data[member_indices[min_idx]] # 原始空间中的代表点
centers.append(center_point)
centers = np.array(centers)
# 可视化结果
plt.figure(figsize=(8, 5))
for j in range(3):
cluster = data[labels == j]
plt.scatter(cluster[:, 0], cluster[:, 1], label=f'聚类{j + 1}')
# 绘制聚类中心
if len(centers) > 0:
plt.scatter(centers[:, 0], centers[:, 1], marker='x', s=100, c='black',
linewidths=2, label='聚类中心', zorder=5)
plt.title("核K均值聚类结果(RBF核 γ=0.5)")
plt.xlabel("花萼长度")
plt.ylabel("花瓣长度")
plt.legend()
plt.tight_layout()
plt.show()
```
### 修改说明:
1. **返回值修改**:`kernel_kmeans` 函数现在返回核矩阵 $K$,用于后续计算中心点。
2. **中心点计算**:
- 遍历每个簇,通过核矩阵计算样本到特征空间中心的距离。
- 选择距离最小的样本作为原始空间中的代表点。
3. **可视化增强**:使用黑色"x"标记中心点,`zorder=5` 确保标记显示在最上层。
### 关键公式说明:
在特征空间中,样本 $x_i$ 到簇中心 $\mu_j$ 的平方距离为:
$$ \|\phi(x_i) - \mu_j\|^2 = K(x_i, x_i) - \frac{2}{m_j}\sum_{k \in C_j} K(x_i, x_k) + \frac{1}{m_j^2}\sum_{k,l \in C_j} K(x_k, x_l) $$
其中 $C_j$ 表示第 $j$ 个簇,$m_j$ 为簇内样本数。通过最小化该距离,找到代表点。
此修改通过核技巧间接获得可视化的聚类中心,增强了结果解释性。
阅读全文
相关推荐














