相關系數



title: 相關系數
date: 2020-01-27 11:42:46
categories: 數學建模
tags: [統計, MATLAB, spss]
mathjax: true

學習視頻:【強烈推薦】清風:數學建模算法、編程和寫作培訓的視頻課程以及Matlab

老師講得很詳細,很受用!!!

相關系數(皮爾遜相關系數)

(1)如果兩個變量本身就是線性的關系,那么皮爾遜相關系數絕對值大的就是相關性強,小的就是相關性弱;
(2)在不確定兩個變量是什么關系的情況下,即使算出皮爾遜相關系數,發現很大,也不能說明那兩個變量線性相關,甚至不能說他們相關,我們一定要畫出散點圖(spss很方便)來看才行。

事實上,比起相關系數的大小,我們往往更關注的是顯著性。(假設檢驗)

假設檢驗

對皮爾遜相關系數進行假設檢驗

  1. 第一步:提出原假設H0和備擇假設H1

    兩個假設是截然相反,假設我們計算出了一個皮爾遜相關系數r,我們想檢驗它是否顯著的異於0.那么我們可以這樣設定原假設和備擇假設:H0:r=0,H1:r≠0

  2. 第二步:在原假設成立的條件下,利用我們要檢驗的量構造出一個符合某一分布的統計量

    (注1:統計量相當於我們要檢驗的量的一個函數,里面不能有其他的隨機變量)
    (注2:這里的分布一般有四種:標准正態分布、t分布、x2分布和F分布)

    對於皮爾遜相關系數r而言,在滿足一定條件下,我們可以構造統計量:

    \[t = r\sqrt{(n-2)/(1-r^{2})} \]

    可以證明t是服從自由度為n-2的t分布

  3. 第三步:將我們要檢驗的這個值帶入這個統計量中,可以得到一個特定的值(檢驗值)。
    假設我們現在計算出來的相關系數為0.5,樣本為30,那么我們帶入第二步公式可以得到\(t^{*}=3.05505\)

  4. 第四步:由於我們知道統計量的分布情況,因此我們可以畫出該分布的概率密度函數pdf,並給定一個置信水平,根據這個置信水平查表找到臨界值,並畫出檢驗統計量的接受域和拒絕域。
    例如,我們知道上述統計量服從自由度為28的t分布,代碼

    x = -4:0.1:4;
    y = tpdf(x,28);
    plot(x,y,'-')
    grid on % 在畫出的圖上加上網格線
    

    常見的置信水平有三個:90%,95%和99%,其中95%是三者中最為常用的。
    因為我們這里是雙側檢驗,所以我們需要找出能覆蓋0.95概率的部分t分布表:https://wenku.baidu.com/view/d94dbd116bd97f192279e94a.html

    查表可知,對應的臨界值為2.048,因此我們可以做出接受域和拒絕域。

  5. 第五步:看我們計算出來的檢驗值是落在了拒絕域還是接受域,並下結論。因為我們得到的t*=3.05505>2.048,因此我們可以下結論:
    在95%的置信水平上,我們拒絕原假設H0:r=0,因此r是顯著的不為0的。

p值判斷法

我們得到的檢驗值t*=3.05505,根據這個值,我們可以計算出其對應的那個概率。

disp('該檢驗值對應的p值為: ')
disp((1-tcdf(3.055,28))*2)
%雙側檢驗的p值要乘以2
  • p<0.01,說明在99%的置信水平上拒絕原假設;p值>0.01,說明在99%的置信水平無法拒絕原假設;
  • p<0.05,說明在95%的置信水平上拒絕原假設;p>0.05,說明在95%的置信水平上無法拒絕原假設:
  • p<0.10,說明在90%的置信水平上拒絕原假設;p>0.10,說明在90%的置信水平上無法拒絕原假設;

小補充: \(0.5、 0.5* 、 0.5**、 0.5***\)的含義是什么?(顯著性標記)

  • 0.5*:90%的水平上顯著異於0,p<0.1
  • 0.5**:95%的水平上顯著異於0,p<0.05
  • 0.5***:99%的水平上顯著異於0,p<0.01
%% 計算各列之間的相關系數以及p值
[R,P] = corrcoef(Test)
% 在EXCEL表格中給數據右上角標上顯著性符號吧
P < 0.01 % 標記3顆星的位置
(P < 0.05) .* (P > 0.01) % 標記2顆星的位置
(P < 0.1) .* (P > 0.05) % 標記1顆星的位置

皮爾遜相關系數假設檢驗的條件

第一, 實驗數據通常假設是成對的來自於正態分布的總體(需要驗證)。 因為我們在求皮爾遜相關性系數以后,通常還會用t檢驗之類的方法來進行皮爾遜相關性系數檢驗,而t檢驗是基於數據呈正態分布的假設的。
第二, 實驗數據之間的差距不能太大。 皮爾遜相關性系數受異常值的影響比較
大。
第三:每組樣本之間是獨立抽樣的。 構造t統計量時需要用到。

正態分布JB檢驗(大樣本 n>30)

對於一個隨機變量{Xi},假設其偏度為s,峰度為K,那么我們可以構造JB統計量:

\[JB=\frac{n}{6}[S^{2}+\frac{(k-3)^{2} }{4}] \]

可以證明,如果{Xi}是正態分布,那么在大樣本情況下JB-x(2)(自由度為2的卡方分布)

正態分布的偏度為0,峰度為3

那么進行假設檢驗的步驟如下:
H0:該隨機變量服從正態分布

H1:該隨機變量不服從正態分布

然后計算該變量的偏度和峰度,得到檢驗值JB*,並計算出其對應的p值

將p值與0.05比較,如果小於0.05則可拒絕原假設,否則我們不能拒絕原假設。

MATLAB中進行JB檢驗的語法: [h,p] = jbtest(x,alpha)
當輸出h等於1時,表示拒絕原假設; h等於0則代表不能拒絕原假設。
alpha就是顯著性水平,一般取0.05,此時置信水平為1‐0.05=0.95
x就是我們要檢驗的隨機變量,注意這里的x只能是向量。

%% 假設檢驗部分
x = -4:0.1:4;
y = tpdf(x,28);  %求t分布的概率密度值 28是自由度  
figure(1)
plot(x,y,'-')
grid on  % 在畫出的圖上加上網格線
hold on  % 保留原來的圖,以便繼續在上面操作
% matlab可以求出臨界值,函數如下
tinv(0.975,28)    %    2.0484
% 這個函數是累積密度函數cdf的反函數
plot([-2.048,-2.048],[0,tpdf(-2.048,28)],'r-')
plot([2.048,2.048],[0,tpdf(2.048,28)],'r-')


%% 計算p值
x = -4:0.1:4;
y = tpdf(x,28);
figure(2)
plot(x,y,'-')
grid on 
hold on
% 畫線段的方法
plot([-3.055,-3.055],[0,tpdf(-3.055,28)],'r-')
plot([3.055,3.055],[0,tpdf(3.055,28)],'r-')
disp('該檢驗值對應的p值為:')
disp((1-tcdf(3.055,28))*2)  %雙側檢驗的p值要乘以2

%% 計算各列之間的相關系數以及p值
[R,P] = corrcoef(Test)
% 在EXCEL表格中給數據右上角標上顯著性符號吧
P < 0.01  % 標記3顆星的位置
(P < 0.05) .* (P > 0.01)  % 標記2顆星的位置
(P < 0.1) .* (P > 0.05) % % 標記1顆星的位置
% 也可以使用Spss操作哦 看我演示

%% 正態分布檢驗
% 正態分布的偏度和峰度
x = normrnd(2,3,100,1);   % 生成100*1的隨機向量,每個元素是均值為2,標准差為3的正態分布
skewness(x)  %偏度
kurtosis(x)  %峰度
qqplot(x)
    
% 檢驗第一列數據是否為正態分布
[h,p] = jbtest(Test(:,1),0.05)%只能一組一組檢驗
[h,p] = jbtest(Test(:,1),0.01)

% 用循環檢驗所有列的數據
n_c = size(Test,2);  % number of column 數據的列數
H = zeros(1,6);  % 初始化節省時間和消耗
P = zeros(1,6);
for i = 1:n_c
    [h,p] = jbtest(Test(:,i),0.05);
    H(i)=h;
    P(i)=p;
end
disp(H)
disp(P)


  ### 小樣本3≤n≤50:Shapiro-wilk檢驗  

H0:該隨機變量服從正態分布

H1:該隨機變量不服從正態分布

計算出威爾克統計量后,得到相應的p值

將p值與0.05比較,如果小於0.05則可拒絕原假設,否則我們不能拒絕原假設。

Q-Q圖 (大樣本 )

在統計學中, Q‐Q圖(Q代表分位數Quantile)是一種通過比較兩個概率分布的分位數對這兩個概率分布進行比較的概率圖方法。
首先選定分位數的對應概率區間集合,在此概率區間上,點(x,y)對應於第一個分布的一個分位數x和第二個分布在和x相同概率區間上相同的分位數。
這里,我們選擇正態分布和要檢驗的隨機變量,並對其做出QQ圖,可想而知,如果要檢驗的隨機變量是正態分布,那么QQ圖就是一條直線。
要利用Q‐Q圖鑒別樣本數據是否近似於正態分布,只需看Q‐Q圖上的點是否近似地在一條直線附近。(要求數據量非常大)

斯皮爾曼spearman相關系數

定義:X和Y為兩組數據,其斯皮爾曼(等級)相關系數:

\[r_{s}=1-\frac{6\sum_{i=1}^{n}d^{2}_{i} }{n(n^2-1)} \]

\[d_{i}為X_{i}和Y_{i}之間的等級差。 \]

(一個數的等級,就是將它所在的一列數按照從小到大排序后,這個數所在的位置)可以證明:\(r_{s}\)位於-1和1之間。

X Y X的等級 Y的等級 等級差 等級差的平方
3 5 2 1 1 1
8 10 5 4.5 0.5 0.25
4 8 3 3 0 0
7 10 4 4.5 ‐0.5 0.25
2 6 1 2 ‐1 1

注:如果有的數值相同,則將它們所在的位置取算術平均

斯皮爾曼相關系數被定義成等級之間的皮爾遜相關系數。

MATLAB中計算

兩種用法
(1) corr(X , Y , 'type' , 'Spearman')
這里的X和Y必須是列向量哦~
(2) corr(X , 'type' , 'Spearman')
這時計算X矩陣各列之間的斯皮爾曼相關系數

%% 斯皮爾曼相關系數
X = [3 8 4 7 2]'  % 一定要是列向量哦,一撇'表示求轉置
Y = [5 10 9 10 6]'
% 第一種計算方法
1-6*(1+0.25+0.25+1)/5/24

% 第二種計算方法
coeff = corr(X , Y , 'type' , 'Spearman')
% 等價於:
RX = [2 5 3 4 1]
RY = [1 4.5 3 4.5 2]
R = corrcoef(RX,RY)

% 計算矩陣各列的斯皮爾曼相關系數
R = corr(Test, 'type' , 'Spearman')

% 大樣本下的假設檢驗
% 計算檢驗值
disp(sqrt(590)*0.0301)
% 計算p值
disp((1-normcdf(0.7311))*2) % normcdf用來計算標准正態分布的累積概率密度函數

% 直接給出相關系數和p值
[R,P]=corr(Test, 'type' , 'Spearman')

斯皮爾曼相關系數的假設檢驗

小樣本

小樣本情況,即n≤30時,直接查臨界值表即可。

樣本相關系數r必須大於等於表中的臨界值,才能得出顯著的結論。

\[H_{0}:r_{s}=0 \]

\[H_{1}:r_{s}!=0 \]

大樣本

大樣本情況下,統計量\(r_{s}\sqrt{n-1}~N(0,1)\)

\[H_{0}:r_{s}=0 \]

\[H_{1}:r_{s}!=0 \]

我們計算檢驗值\(r_{s}\sqrt{n-1}\),並求出對應的p值與0.05相比即可。

p值大於0.05,因此我們無法拒絕原假設。(和0沒有顯著的差異

% 大樣本下的假設檢驗
% 計算檢驗值
disp(sqrt(590)*0.0301)
% 計算p值
disp((1-normcdf(0.7311))*2) % normcdf用來計算標准正態分布的累積概率密度函數
% 直接給出相關系數和p值
[R,P]=corr(Test, 'type' , 'Spearman')

斯皮爾曼相關系數和皮爾遜相關系數選擇:

  1. 連續數據,正態分布,線性關系,用pearson相關系數是最恰當,當然用spearman相關系數也可以, 就是效率沒有pearson相關系數高。

  2. 上述任一條件不滿足,就用spearman相關系數,不能用pearson相關系數。

  3. 兩個定序數據之間也用spearman相關系數,不能用pearson相關系數。

    定序數據是指僅僅反映觀測對象等級、順序關系的數據,是由定序尺度計量形成的,表現為類別,可以進行排序,屬於品質數據。
    例如:優、良、差;
    我們可以用1表示差、 2表示良、 3表示優,但請注意,用2除以1得出的2並不代表任何含義。定序數據最重要的意義代表了一組數據中的某種邏輯順序。
    注:斯皮爾曼相關系數的適用條件比皮爾遜相關系數要廣,只要數據滿足單調關系(例如線性函數、指數函數、對數函數等)就能夠使用

總代碼

clear;clc
load 'physical fitness test.mat'  %文件名如果有空格隔開,那么需要加引號
% https://ww2.mathworks.cn/help/matlab/ref/corrcoef.html
%% 統計描述
MIN = min(Test);  % 每一列的最小值
MAX = max(Test);   % 每一列的最大值
MEAN = mean(Test);  % 每一列的均值
MEDIAN = median(Test);  %每一列的中位數
SKEWNESS = skewness(Test); %每一列的偏度
KURTOSIS = kurtosis(Test);  %每一列的峰度
STD = std(Test);  % 每一列的標准差
RESULT = [MIN;MAX;MEAN;MEDIAN;SKEWNESS;KURTOSIS;STD]  %將這些統計量放到一個矩陣中中表示


%% 計算各列之間的相關系數
% 在計算皮爾遜相關系數之前,一定要做出散點圖來看兩組變量之間是否有線性關系
% 這里使用Spss比較方便: 圖形 - 舊對話框 - 散點圖/點圖 - 矩陣散點圖

R = corrcoef(Test)   % correlation coefficient


%% 假設檢驗部分
x = -4:0.1:4;
y = tpdf(x,28);  %求t分布的概率密度值 28是自由度  
figure(1)
plot(x,y,'-')
grid on  % 在畫出的圖上加上網格線
hold on  % 保留原來的圖,以便繼續在上面操作
% matlab可以求出臨界值,函數如下
tinv(0.975,28)    %    2.0484
% 這個函數是累積密度函數cdf的反函數
plot([-2.048,-2.048],[0,tpdf(-2.048,28)],'r-')
plot([2.048,2.048],[0,tpdf(2.048,28)],'r-')


%% 計算p值
x = -4:0.1:4;
y = tpdf(x,28);
figure(2)
plot(x,y,'-')
grid on 
hold on
% 畫線段的方法
plot([-3.055,-3.055],[0,tpdf(-3.055,28)],'r-')
plot([3.055,3.055],[0,tpdf(3.055,28)],'r-')
disp('該檢驗值對應的p值為:')
disp((1-tcdf(3.055,28))*2)  %雙側檢驗的p值要乘以2

%% 計算各列之間的相關系數以及p值
[R,P] = corrcoef(Test)
% 在EXCEL表格中給數據右上角標上顯著性符號吧
P < 0.01  % 標記3顆星的位置
(P < 0.05) .* (P > 0.01)  % 標記2顆星的位置
(P < 0.1) .* (P > 0.05) % % 標記1顆星的位置
% 也可以使用Spss操作哦 看我演示

%% 正態分布檢驗
% 正態分布的偏度和峰度
x = normrnd(2,3,100,1);   % 生成100*1的隨機向量,每個元素是均值為2,標准差為3的正態分布
skewness(x)  %偏度
kurtosis(x)  %峰度
qqplot(x)
    
% 檢驗第一列數據是否為正態分布
[h,p] = jbtest(Test(:,1),0.05)%只能一組一組檢驗
[h,p] = jbtest(Test(:,1),0.01)

% 用循環檢驗所有列的數據
n_c = size(Test,2);  % number of column 數據的列數
H = zeros(1,6);  % 初始化節省時間和消耗
P = zeros(1,6);
for i = 1:n_c
    [h,p] = jbtest(Test(:,i),0.05);
    H(i)=h;
    P(i)=p;
end
disp(H)
disp(P)

% Q-Q圖
qqplot(Test(:,1))

%% 斯皮爾曼相關系數
X = [3 8 4 7 2]'  % 一定要是列向量哦,一撇'表示求轉置
Y = [5 10 9 10 6]'
% 第一種計算方法
1-6*(1+0.25+0.25+1)/5/24

% 第二種計算方法
coeff = corr(X , Y , 'type' , 'Spearman')
% 等價於:
RX = [2 5 3 4 1]
RY = [1 4.5 3 4.5 2]
R = corrcoef(RX,RY)

% 計算矩陣各列的斯皮爾曼相關系數
R = corr(Test, 'type' , 'Spearman')

% 大樣本下的假設檢驗
% 計算檢驗值
disp(sqrt(590)*0.0301)
% 計算p值
disp((1-normcdf(0.7311))*2) % normcdf用來計算標准正態分布的累積概率密度函數

% 直接給出相關系數和p值
[R,P]=corr(Test, 'type' , 'Spearman')

作業

寫一篇文章,分析男生體測數據各指標之間的相關性,並與女生的數據得到的結論進行對比。
要求:要說明選擇哪一種相關系數的原因,並要求做出散點圖。(可以自己動手試試相關矩陣可視化哦)

確定相關系數

​ 為了度量兩個變量間的線性關系,一般采用皮爾遜(Pearson)相關系數或斯皮爾曼(Spearman)相關系數進行分析。其中,只有當數據滿足連續且呈線性關系時,才能使用皮爾遜(Pearson)相關系數描述變量間的相關性,且估計數據顯著性時候需要數據服從正態分布;相對地,斯皮爾曼(Spearman)相關系數則沒有使用條件限制。
​ 因此,在確定使用哪種相關系數之前,應對各指標數據繪制散點圖判斷是否呈線性關系以及正態分布檢驗。

繪制散點圖

由散點圖司知,各指標彼此之間均無線性關系。

正態分布檢驗

由於本整數據的樣本容里;730,屬於大樣本苷量,應采用JB檢始的方式檢驗備描標款據是否願從正態分布·我們將通過MATLAB 的 jbtest的數對各指標數據進行正態分布檢驗。
假設原指標均服從正態分布,在95%的置信水平下,各指標數據的正態分布檢驗結果如下表所示:
 ![](https://img2018.cnblogs.com/blog/1880713/202002/1880713-20200220234314812-591975900.png)


由表1可知,經正態分布檢驗之后,各指標的h值均為1且p值均小於0.05,即拒絕原假設。
經散點圖與正態分布檢驗分析可知,本不能使用皮爾遜(Pearson)相關系數分析.故考慮使用斯皮爾曼(Speaman)相關系數。

斯皮爾曼(Spearman)相關系數

我們將本題數據導入到SPSS中,進行相關性分析,結果如下表所示:

結果分析。

根據SPSS求解結果,可得結論如下:
當顯著性水平a=0.01時,身高與肺活量呈正相關,身高與坐位體前屈呈負相關,體重與立定跳遠呈正相關。

附錄

%% 導入電子表格中的數據
% 用於從以下電子表格導入數據的腳本:
%
%    工作簿: G:\數學建模學習材料\參考資料\清風數學建模\第1-14講和番外篇的課件和代碼(1月16日修訂版本)\第1-14講和番外篇課件和代碼\第5講.相關系數\代碼和例題數據\八年級男生體測數據.xls
%    工作表: 八年級數據
%
% 要擴展代碼以供其他選定數據或其他電子表格使用,請生成函數來代替腳本。

% 由 MATLAB 自動生成於 2020/02/20 23:38:50

%% 導入數據
[~, ~, raw] = xlsread('G:\數學建模學習材料\參考資料\清風數學建模\第1-14講和番外篇的課件和代碼(1月16日修訂版本)\第1-14講和番外篇課件和代碼\第5講.相關系數\代碼和例題數據\八年級男生體測數據.xls','八年級數據','A2:F731');

%% 創建輸出變量
data = reshape([raw{:}],size(raw));

%% JB檢驗
n_c = size(data,2);  % number of column 數據的列數
H = zeros(1,6);  % 初始化節省時間和消耗
P = zeros(1,6);
for i = 1:n_c
    [h,p] = jbtest(data(:,i),0.05);
    H(i)=h;
    P(i)=p;
end
disp(H)
disp(P)



%% 清除臨時變量
clearvars data raw;


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM