时域上的Pearson相关系数Matlab函数（计算两个变量随时间变化的相关系数）

写在前面

在数据处理中，我们常常需要判断两个时间序列在时域上的相关性，Pearson相关系数广泛用于度量两个变量之间的相关程度，它是由Karl Pearson从Francis Galton在19世纪80年代提出的一个相似却又稍有不同的想法演变而来的。这里首先简单介绍了 Pearson相关系数，接着给出了一个计算两个变量的Pearson相关系数随时间变化的例子，最后给出了一个可以设置时间窗的计算两个变量的Pearson相关系数随时间变化的函数。

Pearson相关系数的定义

两个变量之间的Pearson相关系数定义为两个变量之间的协方差和标准差的商，变量X和Y的Pearson相关系数为：

$${\rho }_{X,Y}=\frac{cov\left(X,Y\right)}{{\sigma }_X{\sigma }_Y}=\frac{E\left[\left(X-{\mu }_X\right)\left(X-{\mu }_Y\right)\right]}{{\sigma }_X{\sigma }_Y}.$$

Pearson相关系数衡量的是线性相关关系。若ρ=0，只能说x与y之间无线性相关关系，不能说无相关关系。相关系数的绝对值越大，相关性越强：相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。
对于x,y之间的相关系数ρ：
当ρ大于0小于1时表示x和y正相关关系
当ρ大于-1小于0时表示x和y负相关关系
当ρ=1时表示x和y完全正相关，ρ=-1表示x和y完全负相关
当ρ=0时表示x和y不相关
通常情况下通过以下取值范围判断变量的相关强度：
相关系数 相关强度
0.8-1.0 极强相关
0.6-0.8 强相关
0.4-0.6 中等程度相关
0.2-0.4 弱相关
0.0-0.2 极弱相关或无相关

链接:
Pearson相关系数-百度百科
皮尔逊相关系数-百度百科

Matlab中的相关函数corrcoef和corr

R = corrcoef(x,y)表示计算序列x和序列y的相关系数，得到的结果是一个2*2矩阵，其中对角线上的元素分别表示x和y的自相关系数，非对角线上的元素分别表示x与y的相关系数和y与x的相关系数。

rho = corr(x,y)表示计算x和y的相关系数矩阵，若x, y是一维列向量，则返回它们的相关系数。
Tipscorr函数有个参数’type’，默认值是’Pearson’

Value	Description
‘Pearson’	Pearson’s Linear Correlation Coefficient
‘Kendall’	Kendall’s Tau Coefficient
‘Spearman’	Spearman’s Rho

计算两个变量Pearson相关系数随时间变化的例子

这里我们给出一个使用Matlab计算两个变量的相关性随时间变化的例子。
定义两个变量x和y，时间长度为1天，采样率为10Hz，则每个变量有24×60×60×10=864000个点。
我们令x为-1到1的随机变量。y分为四个部分：

1. 第一部分（1:21600）与对应x的值相同；
2. 第二部分（21601:432000）与对应x的值正负相反；
3. 第三部分（432001:648000）-1到1的随机变量；
4. 第四部分（648001:864000）是对应x的翻转（顺序颠倒）。

假设以10秒为时间窗的长度来计算一次Pearson相关系数，当我们将时间分段计算时，其实有两种方式：
第一种方式，每个时间窗之间不重叠，则一共要计算86400/10=8640个Pearson相关系数。
第二种方式，每个时间窗之间有重叠，若重叠50%，则除了最开始的5秒和最后的5秒，每个时间窗都被重复计算了一次，一共要计算(86400/5)-1=17279个Pearson相关系数。

这里我们采用每10秒计算一个x与y的Pearson相关系数，相邻时间段之间重叠一半。
结果如下图：

代码如下：

% A example of calculating the Pearson's Linear Correlation Coefficient of 
% two variables x and y in time domain
% by Yuechu Wu, 2021-05-28

clear all;

FS = 10;
npoints = 24*60*60*FS;
part = npoints/4;
pearcorr=[];

% Define variables x and y
x = -1+(1-(-1))*rand(npoints,1);
y = zeros(npoints,1); % Initialization variable y
y(1:part) = x(1:part);
y(part+1:2*part) = -x(part+1:2*part);
y(2*part+1:3*part) = -1+(1-(-1))*rand(part,1);
y(3*part+1:4*part) = flipud(x(3*part+1:4*part));

timeWindow = 10*FS;
step = timeWindow/2;
itera = (npoints/step)-1;

for i = 1:itera
    xi = x((i-1)*step+1:(i+1)*step);
    yi = y((i-1)*step+1:(i+1)*step);
    peari=corr(xi,yi);
    pearcorr = [pearcorr;peari];
end

pearend = corr(x(),y());
pearcorr = [pearcorr;pearend];


figure(1);

subplot(3,1,1)
plot(x(1:step:end),'color',[0,0,139]/255,'Linewidth',1);
xlim([0 17280]);ylim([-1 1]);
xticks(0:1440:17280);
xticklabels({'00','02','04','06','08','10','12','14','16','18','20','22','24'});
xlabel('Time (hour)','FontSize',12);
ylabel('Value of x','FontSize',12);

subplot(3,1,2)
plot(y(1:step:end),'color',[139,0,139]/255,'Linewidth',1);
xlim([0 17280]);ylim([-1 1]);
xticks(0:1440:17280);
xticklabels({'00','02','04','06','08','10','12','14','16','18','20','22','24'});
xlabel('Time (hour)','FontSize',12);
ylabel('Value of y','FontSize',12);

subplot(3,1,3)
plot(pearcorr,'k','Linewidth',2);
xlim([0 17280]);ylim([-1 1]);
xticks(0:1440:17280);
xticklabels({'00','02','04','06','08','10','12','14','16','18','20','22','24'});
xlabel('Time (hour)','FontSize',12);
ylabel('Pearson Correlation');

计算时域上的Pearson相关系数的函数

为了能够能方便的计算两个变量在时域上的相关性，这里给出可以设置时间窗的计算两个变量的Pearson相关系数随时间变化的函数。
函数如下：

function PearCorr = CalcPearson(x,y,timeWindow)

% calculating the Pearson's Linear Correlation Coefficient of 
% two variables x and y in time domain
% x and y are one-dimensional column vectors of the same length
% timeWindow is the number of points to calculate a correlation coefficient
% Yuechu Wu
% 2021-05-28

PearCorr = [];
npoints = length(x);
step = timeWindow/2;
itera = (npoints/step)-1;

for i = 1:itera
    xi = x((i-1)*step+1:(i+1)*step);
    yi = y((i-1)*step+1:(i+1)*step);
    peari=corr(xi,yi);
    PearCorr = [PearCorr;peari];
end

pearend = corr(x(end-step+1:end),y(end-step+1:end));
PearCorr = [PearCorr;pearend];

end

这时，我们回到刚才的例子，使用20秒作为时间窗，并且通过调用刚刚定义的CalcPearson函数来实现计算x和y的相关性，就显得方便很多。
结果如下图：

代码如下：

% A example of calculating the Pearson's Linear Correlation Coefficient of 
% two variables x and y in time domain
% by Yuechu Wu, 2021-05-28

clear all;

FS = 10;
npoints = 24*60*60*FS;
part = npoints/4;

%Define variables x and y
x = -1+(1-(-1))*rand(npoints,1);
y = zeros(npoints,1); %Initialization variable y
y(1:part) = x(1:part);
y(part+1:2*part) = -x(part+1:2*part);
y(2*part+1:3*part) = -1+(1-(-1))*rand(part,1);
y(3*part+1:4*part) = flipud(x(3*part+1:4*part));

timeWindow = 20*FS;
step = timeWindow/2;
pearcorr = CalcPearson(x,y,timeWindow);

figure(1);

subplot(3,1,1)
plot(x(1:step:end),'color',[0,0,139]/255,'Linewidth',1);
ylabel('Value of x');

subplot(3,1,2)
plot(y(1:step:end),'color',[139,0,139]/255,'Linewidth',1);
ylabel('Value of y');

subplot(3,1,3)
plot(pearcorr,'k','Linewidth',2);
ylabel('Pearson Correlation');

问题

• 问题一
  时间窗设置为偶数，这就导致它没有中间数。例如在第一个例子中，第一个时间窗为第1到第100个点，第二个时间窗为第51到第150个点，这就导致了两个时间窗的overlap并不正好是50%。
• 问题二
  当我们计算完时域上的相关系数之后，常常需要将两个变量曲线和相关性曲线画在同一个坐标系上。变量的长度显然要大于相关系数的长度，相关系数的长度为(npoints/step)-1，这个-1有点讨厌，因为我们希望它们是一个整数倍数的关系。当然也可以通过对相关系数进行插值来解决这个问题，但是这里我们采用了另一个思路。
  我们计算相关系数的时候，最后加了一个点，把最后一个时间窗没有重叠的那部分又计算了一个相关系数，这样我们只需要简单对变量x和y降采样就可以将它们画在同一个坐标系上了。

参考

Pearson相关系数-百度百科
皮尔逊相关系数-百度百科
[1]: https://baike.baidu.com/item/Pearson相关系数/6243913
[2]: https://baike.baidu.com/item/皮尔逊相关系数/12712835