協方差/相關矩陣/相關系數

通過兩組統計數據計算而得的協方差可以評估這兩組統計數據的相似程度。

樣本：

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()