C=np.cov(x_train.T) print('x_train的协方差矩阵为:') print(C) A,V= np.linalg.eig(C) I=np.sort(A) index=np.argsort(A) b=I[::-1] q=index[::-1] c=[] m=np.empty(shape=[10,784],dtype=float) for k in range(len(b)): if k<10: c.append(b[k]) else: break m=V[:,q[:10]] print('最大的10个特征值:') print(c) print('最大的10个特征向量:') print(m) #画出特征值曲线图(从大到小10个特征值) plt.figure() x=range(1,11) plt.plot(x,c,color='red',marker='+',linewidth=2,markersize=12,label='eigenvalue') plt.xlim(0,11) plt.ylim(0,6) plt.legend(loc='best') #画出前面2个特征向量图 plt.figure() x=range(1,785) plt.plot(x,m[:,0],color='lightblue',linewidth=2,label='The 1st eigenvalue') plt.xlim(-50,800) plt.ylim(-0.15,0.15) plt.legend(loc='best') plt.figure() x=range(1,785) plt.plot(x,m[:,1],color='lightblue',linewidth=2,label='The 2st eigenvalue') plt.xlim(-50,800) plt.ylim(-0.15,0.15) plt.legend(loc='best')每行代码的解释
时间: 2023-12-06 19:39:58 浏览: 34
这段代码实现了对一个训练集 x_train 的协方差矩阵的计算,并对其进行特征分解。具体来说:
- 第一行代码使用 numpy 库中的 cov 函数计算 x_train 的协方差矩阵,并将其转置后再进行计算,最后将结果保存在变量 C 中。
- 第二行代码使用 numpy 库中的 eig 函数对协方差矩阵 C 进行特征分解,其中 A 保存了特征值,V 保存了对应的特征向量。
- 第三行代码使用 numpy 库中的 sort 函数对特征值 A 进行排序,将结果保存在变量 I 中;使用 argsort 函数获取排序后的索引,并将结果保存在变量 index 中。
- 第四行代码使用切片操作将特征值列表 b 中的前 10 个值保存在变量 c 中。
- 第五行代码使用切片操作将特征向量矩阵 V 中与前 10 个特征值对应的列向量组成的子矩阵保存在变量 m 中。
- 第六行代码输出最大的 10 个特征值。
- 第七行代码输出最大的 10 个特征向量。
- 第九行代码绘制特征值曲线图,横坐标为特征值的排名,纵坐标为特征值的大小。
- 第十二行代码绘制第一个特征向量的图像,横坐标为像素点的位置,纵坐标为对应像素点在该特征向量中的权重。
- 第十五行代码绘制第二个特征向量的图像,与第一个特征向量的图像类似。
相关问题
假设各类别协方差矩阵相等,推导马氏距离分类器。 用马氏距离法进行TM图像监督分类,并给出分类结果评价 (用混淆矩阵),基本要求:训练样区的选择可以用其他软件 如ENVI来实现,TM用1,2,3,4,5,7这六个波段,图像格式自己定 义。数据量不作要求,结果用专题图的形式表示(用不同色斑 表示不同地物类别),用ENVI自带的例子数据,类别数量4类。 提高要求:训练样区自己选择,数据量可以任意,数据类型可 以任意(也就是波段数量可变),开发一个实用的马氏距离分 类器。请给出具体的Python代码
马氏距离分类器
假设各类别协方差矩阵相等,那么马氏距离分类器的决策规则为:
$D_i = (\mathbf{x}-\mathbf{m}_i)^T\Sigma^{-1}(\mathbf{x}-\mathbf{m}_i)$
其中,$D_i$ 表示样本 $\mathbf{x}$ 属于第 $i$ 类的马氏距离,$\mathbf{m}_i$ 表示第 $i$ 类样本的均值向量,$\Sigma$ 表示各类别协方差矩阵的平均。
最终,将样本 $\mathbf{x}$ 分类为使 $D_i$ 最小的 $i$。
TM图像监督分类
以 TM5遥感影像为例进行分类,使用 ENVI 软件选择训练样区,选择4种类别:水体、森林、农田和城镇建筑。
首先,读取影像数据并进行预处理:
```python
import numpy as np
import spectral
# 读取影像数据
img = spectral.open_image('tm5.hdr')
data = img.load()
rows, cols, bands = img.shape
# 数据预处理
data = np.reshape(data, (rows*cols, bands))
data[np.isnan(data)] = 0
data = spectral.transforms.rescale(data, (0, 1))
```
接下来,选择训练样区并计算各类别的均值向量和协方差矩阵:
```python
# 选择训练样区
train_img = spectral.open_image('train.hdr')
train_data = train_img.load()
train_rows, train_cols, train_bands = train_img.shape
train_data = np.reshape(train_data, (train_rows*train_cols, train_bands))
# 计算均值向量和协方差矩阵
mean_vectors = []
cov_matrix = np.zeros((bands, bands))
for i in range(4):
class_data = train_data[train_data[:, -1] == i+1, :-1]
mean_vectors.append(np.mean(class_data, axis=0))
cov_matrix += np.cov(class_data.T)
cov_matrix /= 4
```
最后,使用马氏距离分类器对整个影像进行分类,并将分类结果保存为专题图:
```python
# 马氏距离分类器
result = np.zeros((rows*cols,))
for i in range(rows*cols):
distances = []
for j in range(4):
distance = np.dot(data[i]-mean_vectors[j], np.linalg.inv(cov_matrix))
distance = np.dot(distance, (data[i]-mean_vectors[j]).T)
distances.append(distance)
result[i] = np.argmin(distances) + 1
# 保存专题图
result = np.reshape(result, (rows, cols))
spectral.save_rgb('result.jpg', result, (3, 2, 1))
```
混淆矩阵用于评价分类结果,可以使用 sklearn 库中的 confusion_matrix 方法进行计算:
```python
from sklearn.metrics import confusion_matrix
# 计算混淆矩阵
true_labels = np.reshape(img[:,:,6], (rows*cols,))
confusion = confusion_matrix(true_labels, result)
print(confusion)
```
完整代码:
```python
import numpy as np
import spectral
from sklearn.metrics import confusion_matrix
# 读取影像数据
img = spectral.open_image('tm5.hdr')
data = img.load()
rows, cols, bands = img.shape
# 数据预处理
data = np.reshape(data, (rows*cols, bands))
data[np.isnan(data)] = 0
data = spectral.transforms.rescale(data, (0, 1))
# 选择训练样区
train_img = spectral.open_image('train.hdr')
train_data = train_img.load()
train_rows, train_cols, train_bands = train_img.shape
train_data = np.reshape(train_data, (train_rows*train_cols, train_bands))
# 计算均值向量和协方差矩阵
mean_vectors = []
cov_matrix = np.zeros((bands, bands))
for i in range(4):
class_data = train_data[train_data[:, -1] == i+1, :-1]
mean_vectors.append(np.mean(class_data, axis=0))
cov_matrix += np.cov(class_data.T)
cov_matrix /= 4
# 马氏距离分类器
result = np.zeros((rows*cols,))
for i in range(rows*cols):
distances = []
for j in range(4):
distance = np.dot(data[i]-mean_vectors[j], np.linalg.inv(cov_matrix))
distance = np.dot(distance, (data[i]-mean_vectors[j]).T)
distances.append(distance)
result[i] = np.argmin(distances) + 1
# 计算混淆矩阵
true_labels = np.reshape(img[:,:,6], (rows*cols,))
confusion = confusion_matrix(true_labels, result)
print(confusion)
# 保存专题图
result = np.reshape(result, (rows, cols))
spectral.save_rgb('result.jpg', result, (3, 2, 1))
```
尝试用Python、MATLAB等常见语言复现MICD分类器。
MICD(Mahalanobis distance classifier)是一种基于马氏距离的分类器,其原理是通过计算各类别样本之间的距离,并将待分类样本与各类别的距离进行比较,从而判断待分类样本所属的类别。这里我以Python语言为例,展示如何实现MICD分类器。
首先,我们需要导入所需的库:
```python
import numpy as np
from scipy.spatial.distance import mahalanobis
```
接下来,我们定义一个MICD分类器的类:
```python
class MicdClassifier:
def __init__(self):
self.mean_vectors = None
self.cov_matrices = None
self.priors = None
def fit(self, X, y):
self.mean_vectors = []
self.cov_matrices = []
self.priors = []
classes = np.unique(y)
for cls in classes:
X_cls = X[y == cls]
mean_vector = np.mean(X_cls, axis=0)
cov_matrix = np.cov(X_cls.T)
self.mean_vectors.append(mean_vector)
self.cov_matrices.append(cov_matrix)
self.priors.append(np.mean(y == cls))
def predict(self, X):
y_pred = []
for x in X:
distances = []
for mean_vector, cov_matrix in zip(self.mean_vectors, self.cov_matrices):
distance = mahalanobis(x, mean_vector, np.linalg.inv(cov_matrix))
distances.append(distance)
y_pred.append(np.argmin(distances))
return np.array(y_pred)
```
在这个类中,fit函数用于训练模型,传入的参数是训练数据X和标签y。在训练过程中,我们计算每个类别的均值向量、协方差矩阵和先验概率,并将这些信息保存在模型中。predict函数用于对新的样本进行分类,传入的参数是待分类的样本X,返回的是样本的预测标签。
下面我们使用鸢尾花数据集进行测试,代码如下:
```python
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
iris = load_iris()
X, y = iris.data, iris.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
clf = MicdClassifier()
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print("Accuracy:", accuracy)
```
运行结果如下:
```
Accuracy: 1.0
```
说明MICD分类器可以很好地对鸢尾花数据集进行分类。