我将带头冲锋——典型相关分析

---

联系数计加油站获取论文和代码等资源。
邮箱:MathComputerSharing@outlook.com
微信公众号:数计加油站
致谢:清风数学建模

---

组与组的整体相关性

相关系数能衡量一对数据的线性关系。然而，我们常常会碰到这样的情况：有两组数据，每组各有几列数据，需要我们研究研究这两组数据之间有没有什么难以告人的关系。

一种做法是两两计算两组数据的每两列数据之间的相关系数，这样可以得出两组数据各列之间的关系。但我们想要得到的，往往并不是这种个体与个体之间的关系，而是两大组数据的整体相关性。比如，我们想要知道文科成绩和理科成绩的总体关系，而不仅仅只是语文和数学之间的关系、英语和物理之间的关系等等。

“代表大会制度”的引入

衡量整体相关性，很容易想到的一个方法，就是在每组之中分别找出一个代表，我们来研究这两个代表之间的关系就可以了。代表之间没有关系，那当然各组各选出一列数据，它们的关系也不会大到哪里去；代表之间纠葛很深，那当然各组各选出一列数据，它们的关系也不会那么的形同陌路。

代表之间的关系，表示着从整体上而言，两组数据之间的关系。

怎么挑选代表——尽量产生关系

代表当然不能简单地从已知数据列中挑选。我们可以采用线性组合的方式，以不同的权重，综合数据组中的各列已知数据。

还有一个原则：万事万物一定是有它们的内部联系存在的。我们挑选的两个代表，理应有着最大的关系，这样才有研究的意义。

如果代表之间的最大关系能够非常大，那么在特定代表的前提下，两组数据确实能产生一定联系。

如果代表之间的最大关系还是很小，那么这两组数据在很大程度上，关系很小。

一篇简短的小论文

Matlab代码

clc,clear

%% 导入数据
data = readmatrix("data\健康指标.csv");
[m,n] = size(data);
X = data(:,1:3);
Y = data(:,4:6);
p = size(X,2);
q = size(Y,2);

%% 进行正态性检验
% Shapiro-wilk 检验
sw_H = ones([n,1]);
sw_P = ones([n,1]);
for i=1:n
    [sw_h,sw_p] = swtest(data(:,i),0.1);
    sw_H(i,1)=sw_h;
    sw_P(i,1)=sw_p;
end
disp("    是否通过      p")
disp([sw_H,sw_P]);

%% 未进行标准化的典型相关分析
[A,B,r,U,V,stats] = canoncorr(X,Y);
disp("r=")
disp(r);
disp("p=");
disp(stats.p);
disp("A=")
disp(A);
disp("B=")
disp(B);

%% 进行标准化的典型相关分析
R11 = corr(X,X);
R12 = corr(X,Y);
R21 = corr(Y,X);
R22 = corr(Y,Y);
M = R11\R12/(R22)*R21;
N = R22\R21/(R11)*R12;
[A_normalized,lambda1] = eig(M);
[B_normalized,lambda2] = eig(N);

lambda1 = sum(sqrt(lambda1));
[lambda1,I1] = sort(lambda1,"descend");
A_normalized = A_normalized(:,I1);
lambda2 = sum(sqrt(lambda2));
[lambda2,I2] = sort(lambda2,"descend");
B_normalized = B_normalized(:,I2);

for i=1:p
    d = A_normalized(:,i)' * R11 *A_normalized(:,i);
    A_normalized(:,i)=A_normalized(:,i)/(d^(1/2));
end

for i=1:q
    d = B_normalized(:,i)' * R22 *B_normalized(:,i);
    B_normalized(:,i)=B_normalized(:,i)/(d^(1/2));
end

disp("lambda_1=");
disp(lambda1);
disp("A_normalzed=");
disp(A_normalized);
disp("lambda_2=");
disp(lambda2);
disp("B_normalzed=");
disp(B_normalized);

%% 典型载荷分析
% 未标准化
RUX = corr(U,X);
RUY = corr(U,V);
RVX = corr(V,X);
RVY = corr(V,Y);
% 标准化
nRUX = A_normalized'*R11;
nRUY = A_normalized'*R12;
nRVX = B_normalized'*R21;
nRVY = B_normalized'*R22;

%% 典型冗余分析
RdU = sum(nRUX.^2,2)./p;
RdV = sum(nRVY.^2,2)./q;

clear d I1 I2 i   

Python代码

## 导库
import numpy as np
import pandas as pd
from scipy import stats

np.set_printoptions(suppress=True)

## 数据初始化
data = pd.read_csv("../data/健康指标.csv", header=None)
data.columns = ["体重", "腰围", "脉搏", "引体向上次数", "起坐次数", "跳跃次数"]
display(data)

m, n = data.shape
D = np.array(data)
X = np.array(data[["体重", "腰围", "脉搏"]])
Y = np.array(data[["引体向上次数", "起坐次数", "跳跃次数"]])
p = X.shape[1]
q = Y.shape[1]

## 正态性检验
swtest_result = []
for i in range(n):
    W, sw_p = stats.shapiro(D[:, i])
    swtest_result.append([W, sw_p])
swtest_result = pd.DataFrame(swtest_result, columns=["W值", "p值"], index=data.columns)
display(swtest_result)

def cov(X, Y,normalized=False):
"""计算两个矩阵的协方差矩阵"""
    p = X.shape[1]
    q = Y.shape[1]
    result = np.ones((p, q))
    for i in range(p):
        for j in range(q):
            if normalized:
                result[i, j] = np.corrcoef(X[:, i], Y[:, j])[0, 1]
            else:
                result[i,j] = np.cov(X[:, i], Y[:, j])[0, 1]
    return result

def sortVbyL(L,V):
"""按特征值大小排序特征向量"""
    vectors = V.T.tolist()
    tmp = zip(L,vectors)
    tmp = sorted(tmp,key=lambda x:x[0],reverse=True)
    L,V = zip(*tmp)
    L = np.array(L)
    V = np.array(V).T
    return L,V

def canoncorr(X, Y,normalized=False):
"""典型相关分析"""
    # 计算协方差矩阵
    S11 = cov(X, X,normalized)
    S12 = cov(X, Y,normalized)
    S21 = cov(Y, X,normalized)
    S22 = cov(Y, Y,normalized)
    # 求特征值和特征向量并排序
    M = np.linalg.inv(S11) @ S12 @ np.linalg.inv(S22) @ S21
    N = np.linalg.inv(S22) @ S21 @ np.linalg.inv(S11) @ S12

    lambda_1, A = np.linalg.eig(M)
    lambda_1 = np.sqrt(lambda_1)
    lambda_1, A = sortVbyL(lambda_1,A)

    lambda_2, B = np.linalg.eig(N)
    lambda_2 = np.sqrt(lambda_2)
    lambda_2, B = sortVbyL(lambda_2,B)

    if np.linalg.norm(lambda_1 - lambda_2) > 0.0001:
        return None
    r = lambda_1

    # 限制方差为1
    for i in range(p):
        d1 = A[:, i].T @ S11 @ A[:, i]
        A[:, i] = A[:, i] * (d1 ** (-1 / 2))
    for j in range(q):
        d2 = B[:, j].T @ S22 @ B[:, j]
        B[:, j] = B[:, j] * (d2 ** (-1 / 2))

    # 计算典型向量
    U = X @ A
    V = Y @ B

    # 假设检验
    P = []
    for k in range(r.shape[0]):
        Lambda_k = np.prod([1 - x ** 2 for x in r[k:]])
        m_k = (m - k - 1) - 1 / 2 * (p + q + 1) + np.sum([x ** (-2) for x in r[:k]])
        f_k = (p - k) * (q - k)
        P.append(stats.chi2.sf(-m_k * np.log(Lambda_k), f_k))

    return A, B, r, U, V, P

## 未进行标准化的典型相关分析
A, B, r, U, V, P = canoncorr(X, Y)
display(r, P, A, B)

## 进行标准化的典型相关分析
nA, nB, nr, nU, nV, nP = canoncorr(X, Y,normalized=True)
display(nr, nP, nA, nB)

## 典型载荷分析
# 未标准化
RUX = cov(U,X,normalized=True)
RUY = cov(U,V,normalized=True)
RVX = cov(V,X,normalized=True)
RVY = cov(V,Y,normalized=True)
# 标准化
nRUX = nA.T @ cov(X,X,normalized=True)
nRUY = nA.T @ cov(X,Y,normalized=True)
nRVX = nB.T @ cov(Y,X,normalized=True)
nRVY = nB.T @ cov(Y,Y,normalized=True)
display(RUX)
display(RVY)
display(nRUX)
display(nRVY)

## 典型冗余分析
RdU = np.sum(nRUX**2,1)/p
RdV = np.sum(nRVY**2,1)/q
display(RdU)
display(RdV)