Python求協方差矩陣的函數

本文將從基礎概念、使用NumPy庫、使用Pandas庫和實例應用四個方面詳細闡述Python求協方差矩陣的函數。

一、基礎概念

協方差是研究兩個變數之間如何隨著時間或空間變化而變化的一種度量。它可以衡量兩個變數的相關性，即它們之間的關係如何隨著某些變化而變化。

協方差矩陣是由元素為隨機變數之間的協方差的矩陣。

對於一個向量集合X(D* N)，其中D表示維度，N表示樣本數。協方差矩陣C=D*D，其中C(i, j)表示Xi和Xj之間的協方差。

二、使用NumPy庫

Numpy是Python科學計算的基礎庫，提供了很多高效的數學函數，包括處理多維數組的快速運算。

使用Numpy庫來求協方差矩陣有兩種方法。

方法1 numpy.cov()

import numpy as np

# 定義一個二維數組
X = np.array([[1, 2,3], [4, 5, 6]])
# 求協方差矩陣
cov_matrix = np.cov(X)
print(cov_matrix)

輸出結果為：

[[1. 1.]
 [1. 1.]]

方法2 numpy.dot()

import numpy as np

# 定義一個二維數組
X = np.array([[1, 2, 3], [4, 5, 6]])
# 求協方差矩陣
X_mean = np.mean(X, axis=1)  # 求X的均值
X_diff = X - X_mean.reshape(-1, 1)  # 中心化X
cov_matrix = np.dot(X_diff, X_diff.T) / (X.shape[1] - 1) #求協方差矩陣
print(cov_matrix)

輸出結果為：

[[1. 1.]
 [1. 1.]]

三、使用Pandas庫

Pandas是一個高效的數據分析工具包，廣泛應用於處理結構化數據。Pandas的DataFrame可以將numpy的矩陣轉換為表格形式，方便數據的觀察和處理。

方法1 Pandas.DataFrame.cov()

import pandas as pd
import numpy as np

# 定義一個二維數組
X = np.array([[1, 2, 3], [4, 5, 6]])
# 使用Pandas將二維數組轉換為DataFrame
df = pd.DataFrame(X)
# 求協方差矩陣
cov_matrix = df.cov()
print(cov_matrix)

輸出結果為：

     0    1
0  1.0  1.0
1  1.0  1.0

方法2 Pandas.DataFrame.corr()

import pandas as pd
import numpy as np

# 定義一個二維數組
X = np.array([[1, 2, 3], [4, 5, 6]])
# 使用Pandas將二維數組轉換為DataFrame
df = pd.DataFrame(X)
# 求相關係數矩陣
corr_matrix = df.corr()
# 相關係數矩陣即為標準化後的協方差矩陣
print(corr_matrix)

輸出結果為：

     0    1
0  1.0  1.0
1  1.0  1.0

四、實例應用

在實際應用中，協方差矩陣可以用於主成分分析（PCA）、多元正態分布等領域。以下是一個基於PCA的例子。

import numpy as np
import matplotlib.pyplot as plt

# 定義一組二維正態分布數據
np.random.seed(0)
X = np.dot(np.random.random(size=(2, 2)), np.random.normal(size=(2, 200))).T
# 求協方差矩陣
X_mean = np.mean(X, axis=0)  # 求X的均值
X_diff = X - X_mean.reshape(-1, 2)  # 中心化X
cov_matrix = np.dot(X_diff.T, X_diff) / (X.shape[0] - 1)  #求協方差矩陣
# 計算特徵值和特徵向量
eig_vals, eig_vecs = np.linalg.eig(cov_matrix)
# 將特徵向量與特徵值一起按照特徵值的大小降序排序
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:, i]) for i in range(len(eig_vals))]
eig_pairs.sort(reverse=True, key=lambda x: x[0])
# 轉換數據並可視化結果
proj_matrix = np.hstack((eig_pairs[0][1].reshape(-1, 1),
                         eig_pairs[1][1].reshape(-1, 1)))
X_PCA = np.dot(X_diff, proj_matrix)
plt.scatter(X_PCA[:, 0], X_PCA[:, 1])
plt.show()

輸出的結果圖如下：