协方差/相关矩阵/相关系数

Posted 2022-10-03 maplethefox

通过两组统计数据计算而得的协方差可以评估这两组统计数据的相似程度。

样本：

A = [a1, a2, ..., an]
B = [b1, b2, ..., bn]

平均值：

ave_a = (a1 + a2 +...+ an)/n
ave_b = (b1 + b2 +...+ bn)/m

离差（用样本中的每一个元素减去平均数，求得数据的误差程度）：

dev_a = [a1, a2, ..., an] - ave_a
dev_b = [b1, b2, ..., bn] - ave_b

协方差

协方差可以简单反映两组统计样本的相关性，值为正，则为正相关；值为负，则为负相关，绝对值越大相关性越强。

cov_ab = ave(dev_a x dev_b)
cov_ba = ave(dev_b x dev_a)

案例：计算两组数据的协方差，并绘图观察。

import numpy as np
import matplotlib.pyplot as mp

a = np.random.randint(1, 30, 10)
b = np.random.randint(1, 30, 10)
#平均值
ave_a = np.mean(a)
ave_b = np.mean(b)
#离差
dev_a = a - ave_a
dev_b = b - ave_b
#协方差
cov_ab = np.mean(dev_a*dev_b)
cov_ba = np.mean(dev_b*dev_a)
print(‘a与b数组：‘, a, b)
print(‘a与b样本方差：‘, np.sum(dev_a**2)/(len(dev_a)-1), np.sum(dev_b**2)/(len(dev_b)-1))
print(‘a与b协方差：‘,cov_ab, cov_ba)
#绘图，查看两条图线的相关性
mp.figure(‘COV LINES‘, facecolor=‘lightgray‘)
mp.title(‘COV LINES‘, fontsize=16)
mp.xlabel(‘x‘, fontsize=14)
mp.ylabel(‘y‘, fontsize=14)
x = np.arange(0, 10)
#a,b两条线
mp.plot(x, a, color=‘dodgerblue‘, label=‘Line1‘)
mp.plot(x, b, color=‘limegreen‘, label=‘Line2‘)
#a,b两条线的平均线
mp.plot([0, 9], [ave_a, ave_a], color=‘dodgerblue‘, linestyle=‘--‘, alpha=0.7, linewidth=3)
mp.plot([0, 9], [ave_b, ave_b], color=‘limegreen‘, linestyle=‘--‘, alpha=0.7, linewidth=3)

mp.grid(linestyle=‘--‘, alpha=0.5)
mp.legend()
mp.tight_layout()
mp.show()

相关系数

协方差除去两组统计样本的乘积是一个[-1, 1]之间的数。该结果称为统计样本的相关系数。

# a组样本 与 b组样本做对照后的相关系数
cov_ab/(std_a x std_b)
# b组样本 与 a组样本做对照后的相关系数
cov_ba/(std_b x std_a)
# a样本与a样本作对照   b样本与b样本做对照   二者必然相等
cov_ab/(std_a x std_b)=cov_ba/(std_b x std_a)

通过相关系数可以分析两组数据的相关性：

若相关系数越接近于0，越表示两组样本越不相关。
若相关系数越接近于1，越表示两组样本正相关。
若相关系数越接近于-1，越表示两组样本负相关。

案例：输出案例中两组数据的相关系数。

print(‘相关系数：‘, cov_ab/(np.std(a)*np.std(b)), cov_ba/(np.std(a)*np.std(b)))

相关矩阵

a与a的相关系数　　a与b的相关系数

b与a的相关系数　　b与b的相关系数

矩阵正对角线上的值都为1。（同组样本自己相比绝对正相关）

# 相关矩阵
numpy.corrcoef(a, b)    
# 相关矩阵的分子矩阵 
# [[a方差，ab协方差], [ba协方差, b方差]]
numpy.cov(a, b)        

 

# 协方差
import numpy as np
import matplotlib.pyplot as mp
import datetime as dt
import matplotlib.dates as md


def dmy2ymd(dmy):
  """
  把日月年转年月日
  :param day:
  :return:
  """
  dmy = str(dmy, encoding=‘utf-8‘)
  t = dt.datetime.strptime(dmy, ‘%d-%m-%Y‘)
  s = t.date().strftime(‘%Y-%m-%d‘)
  return s


dates, bhp_closing_prices =   np.loadtxt(‘bhp.csv‘,
             delimiter=‘,‘,
             usecols=(1, 6),
             unpack=True,
             dtype=‘M8[D],f8‘,
             converters=1: dmy2ymd)  # 日月年转年月日
vale_closing_prices =   np.loadtxt(‘vale.csv‘,
             delimiter=‘,‘,
             usecols=(6,),
             unpack=True)  # 因为日期一样,所以此处不读日期
# print(dates)
# 绘制收盘价的折现图
mp.figure(‘COV‘, facecolor=‘lightgray‘)
mp.title(‘COV‘, fontsize=18)
mp.xlabel(‘Date‘, fontsize=14)
mp.ylabel(‘Price‘, fontsize=14)
mp.grid(linestyle=":")

# 设置刻度定位器
# 每周一一个主刻度,一天一个次刻度

ax = mp.gca()
ma_loc = md.WeekdayLocator(byweekday=md.MO)
ax.xaxis.set_major_locator(ma_loc)
ax.xaxis.set_major_formatter(md.DateFormatter(‘%Y-%m-%d‘))
ax.xaxis.set_minor_locator(md.DayLocator())
# 修改dates的dtype为md.datetime.datetiem
dates = dates.astype(md.datetime.datetime)
mp.plot(dates, bhp_closing_prices,
        color=‘dodgerblue‘,
        linewidth=2,
        linestyle=‘-‘,
        alpha=0.8,
        label=‘BHP Closing Price‘)
mp.plot(dates, vale_closing_prices,
        color=‘orangered‘,
        linewidth=2,
        linestyle=‘-‘,
        alpha=0.8,
        label=‘VALE Closing Price‘)
# 计算两只股票收盘价的协方差
# 求均值
mean_bhp = np.mean(bhp_closing_prices)
mean_vale = np.mean(vale_closing_prices)
# 求离差
dev_bhp = bhp_closing_prices - mean_bhp
dev_vale = vale_closing_prices - mean_vale
# 协方差(bhp的离差乘以vale的离差)
cov = (dev_bhp * dev_vale).mean()  # 用的是总体的数据
print(‘cov:‘, cov)  # cov: 3.135577333333333
# 两条图线的相关性
k = cov / np.std(bhp_closing_prices * np.std(vale_closing_prices))
print(k)  # 0.8664988296368301
# 相关系数 = a方差/(a标准差*a标准差)
k = np.corrcoef(bhp_closing_prices, vale_closing_prices)
print(k)
"""
[[1.         0.86649883]
 [0.86649883 1.        ]]
"""
# q秋两组数据的协方差
covm = np.cov(bhp_closing_prices, vale_closing_prices)  # 用的是 样本的数据
print(covm)
"""
[[8.53844379 3.24370069]
 [3.24370069 1.64122023]]
"""

mp.legend()
mp.gcf().autofmt_xdate()
mp.show()