数模第五次作业：自己实现斯皮尔曼相关系数

数模自用

1.斯皮尔曼系数及P值计算

%利用第一种方法计算斯皮尔曼相关系数
%即先对X和Y两个指标的数据进行等级的排序
%并计算等级差
%将其平方
%利用公式计算得出
function [R,P]=fun_spearman(X,kind)

    %判断输入的参数个数，如果为1个，则默认是双侧检验
    %kind=1代表单侧检验，kind=2代表双侧检验
    if nargin==1
        disp("未输入检验种类，默认为双侧检验")
        kind=2;
    end
    
    %分析数据矩阵
    %row：行数 col：列数
    [row,col]=size(X);
    
    %判断样本数
    if row<30
        disp("样本数小于30，请直接查表得出p值")
    elseif col<2
        disp("指标个数小于2，无法计算")
    elseif kind~=1 && kind~=2
        disp("kind的值只能为1或者2，1：单侧检验；2：双侧检验")
    else
        %初始化记录相关系数的矩阵R
        R=zeros(col,col);
        P=zeros(col,col);
        %计算每个指标和另外其他所有指标的相关系数
        for i=1:col
            for j=i+1:col
                %取出每一列和除这一列以外的其他列数据
                x=X(:,i);
                y=X(:,j);
                %计算该指标和另一指标的相关系数
                R(i,j)=calculate_r(x,y);
                %因为是对称矩阵
                R(j,i)=R(i,j);
                P(i,j)=calculate_p(R(i,j),row,kind);
                P(j,i)=P(i,j);
            end
            R(i,i)=1;
            P(i,i)=1;
        end
    end
end

2.计算相关系数

%输入两个指标的列向量，返回相关系数
function [r]=calculate_r(X,Y)

    %分别计算两组数据的rank
    rank_X=calculate_rank(X);
    rank_Y=calculate_rank(Y);
    
    n=size(X,1);

    %计算两组数据的等级差
    D=rank_X-rank_Y;

    %按照第一种方法计算斯皮尔曼相关系数
    r=1-6*sum(D.*D)/(n*(n*n-1));
end

3.计算等级

function [data_rank] = calculate_rank(X)
    %分别计算XY数据的等级
    %先对数据进行排序，得到的index就是按升序排好数据后
    %排好序的数据在原来的矩阵中的索引
    %比如[3 8 4 7 2]排序后变成[2 3 4 7 8]，对应的index为[5 1 3 4 2]
    [~,index]=sort(X);
    %得到排好序后的数据在原矩阵的索引后，还要对索引进行排序，以得到等级
    [~,data_rank]=sort(index);
    
    %得到X数据个数
    row=size(X,1);
    
    %遍历所有数据，确定rank
    for i=1:row
        %返回跟当前数据相等的在矩阵中的数据，相等则sameposition为1
        samePosition=(X==X(i));
        %重新计算data_rank，将相同的数据的rank加起来，除以相同数据的重数
        data_rank(samePosition==1)=sum(data_rank.*samePosition)./sum(samePosition);
    end
end

4.计算p值

%通过相关系数，数据个数，检验的种类计算p值，来判断r是否显著异于0

function [p] = calculate_p(r,n,kind)

    %构造标准正态分布统计量
    z=abs(r)*sqrt(n-1);
    p=(1-normcdf(z))*kind;
end