数学建模（NO.9斯皮尔曼相关系数）

一. 定义

方法一：计算公式法

三.方法二：利用皮尔逊相关系数法

四.Matlab计算

五.斯皮尔曼和皮尔逊对比

六.斯皮尔曼相关系数的假设检验

1.小样本（查表）

2.大样本（计算统计量）

七.两个比较

八.相关代码

%% 斯皮尔曼相关系数
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')

function [p]= calculate_p(r, m, kind)
% % 输入值：
    % r：斯皮尔曼相关系数
    % m: 样本个数
    % kind: 1表示单侧检验 2表示双侧检验
% % 返回值：
    % p：计算出来的p值

    z = abs(r) * sqrt(m-1); % 计算检验值  注意这里的相关系数我们先转换为正数再进行计算
    p = (1 - normcdf(z)) * kind; % 计算p值，双侧检验的p值是单侧检验的2倍

end

function [r]= calculate_r(X, Y)
% % 输入值：
    % X: 列向量
    % Y: 列向量，且与X同维度
% % 返回值：
    % r: X和Y的斯皮尔曼相关系数（第一种定义方法）
    
    RX = rank_data(X);  % 调用自定义函数 rank_data 来计算X的等级
    RY = rank_data(Y);  % 调用自定义函数 rank_data 来计算Y的等级
    d = RX - RY; % 计算X和Y等级差
    n = size(X,1); % 计算样本个数n
    r = 1 - (6 * sum(d .* d)) / (n * (n^2-1));  % 利用公式计算斯皮尔曼相关系数

end

function [ R , P ]= fun_spearman(X, kind)  
% % 输入值：
    % X: m*n维数据矩阵，每一行表示一个样本，每一列表示一个指标；且 m >=30 以及 n >= 2 
    % kind=1: 单侧检验；kind=2: 双侧检验 (不输入默认为2)
% % 返回值：
    % R: 斯皮尔曼相关系数矩阵（n*n维）
    % P: 对应的p值矩阵（n*n维）
    
    if nargin == 1  % 判断用户输入的参数，如果只输入了一个参数，则默认kind = 2
        kind = 2;
    end

    [m,n] = size(X); % 计算样本个数和指标个数

    if m < 30  % 判断是否样本数太少
        disp('样本个数少于30，请直接查临界值表进行假设检验')

    elseif n <2  % 判断是否指标数太少
        disp('指标个数太少，无法计算')

    elseif kind ~= 1 && kind ~= 2 % 判断kind是否为1或者2
        disp('kind只能取1或者2')

    else  % 如果上述输入均没问题的话就执行下面的操作
        R = ones(n); % 初始化R矩阵
        P = ones(n); % 初始化P矩阵
        for i = 1: n
            for j = (i+1): n   % 这样设置循环只计算主对角线上半部分的值
                r = calculate_r(X(:, i), X(:, j));  % 用子函数 calculate_r 计算i和j两列的相关系数r
                p = calculate_p(r, m, kind); % 用子函数 calculate_p 计算p值
                R(i, j) = r;  R(j, i) = r;  % 把计算出来的相关系数r填充到我们的R矩阵中，注意R矩阵对称
                P(i, j) = p; P(j, i) = p;  % 把计算出来的p值填充到我们的P矩阵中，注意P矩阵对称
            end
        end
    end
    
end

function [RX]= rank_data(X)
% % 输入值：
    % X: 列向量
% % 返回值：
    % RX: 对应的X的等级

%  举个例子X = [5 10 9 10 6]'
    [~ ,index] = sort(X);  % ~表示我们不需要第一个输出值（即我们排序后的X [5 6 9 10 10]）
                                      % 注意这里的index = [1 5 3 2 4]' 是我们排序后的X在原向量中的位置
    [~ ,RX] = sort(index);  % 对index进行一次升序，得到的rx就是我们想要的等级 rx = [1 4 3 5 2]'
                                       % 但是这个等级还有一点小问题 ，那就是没有考虑到相等取平均值的问题
                                       
    for i = 1:size(X,1)  % 设置一个循环  (假设此时程序运行到了i = 2)
        position = ( X == X(i) ); % 得到X中与X(i)相等的位置，返回一个列向量，向量值全为1或0
                                              % (i= 2时,position = [0 1 0 1 0]' )
        RX(position == 1) = sum(RX .* position) / sum(position);  % 对RX进行处理 
        %  rx .* position = [0 4 0 5 0]'   
        %  那么 sum(rx .* position) / sum(position) = (4+5) / 2 = 4.5
        %  rx(position == 1)  = 4.5 ： 对rx中与position == 1对应位置的元素（即第2和第4个位置）进行赋值操作
    end
end