matlab滑动t检验和全球海表平均温度突变检验

滑动t检验程序一：

clear all;close all;clc
%%   读取数据
data = xlsread('Tair.xlsx'); % 读取数据
time = data(:,1);%时间序列数据
temperature = data(:,2);%最高气温数据
%%   画时间序列
figure(1)
set(gcf,'color','w')
plot(time,temperature)
xlabel('year')
ylabel('temperature(\circC)')
title('1911-1995年平均气温时间序列')
%%  滑动t检验
step = 10;n1 = 10;n2 = 10;%设置子序列步长
v = step+step-2;%自由度
t=[];
ttest=2.10;%t检验值  alpha=0.05
for i = 1:length(time)-step-step+1
    x1 = temperature(i:i+step-1);
    x2 = temperature(i+step:i+step+step-1);
    meanx1 = mean(x1);
    meanx2 = mean(x2);
    a = meanx1-meanx2;
    b = (n1+n2)/(n1*n2);
    varx1 = var(x1);
    varx2 = var(x2);
    c= ((n1-1)*varx1+(n2-1)*varx2)/(n1+n2-2);
    t1=a/sqrt(c*b);
    t=cat(1,t,t1);
end
x = time(step:length(time)-step);
figure(2)
set(gcf,'color','w')
plot(x,t,'r')
axis([min(x),max(x),min(t),max(t)]);
xlabel('year');
ylabel('t统计量');
hold on
x1=x(:,1);y1=ttest*ones(size(x1));
plot(x1,y1,'b--')
x2=x(:,1);y2=-ttest*ones(size(x2));
plot(x2,y2,'b--')
legend('t统计量','0.05显著水平');
title('1911-1995年平均气温突变检验')

程序二

clear
mydata = xlsread('Tair.xlsx'); % 读取数据
timeSeries = mydata(:,1);%时间序列数据
dataSeries = mydata(:,2);%最高气温数据
dataCount = length(dataSeries);
%%设置步长与检验值
step = 10; %步长
v = step+step-2;  %计算自由度
ttest = 2.878; %查表得t检验值，修改
len1 = step;
len2 = step;
x = timeSeries(step:dataCount-step)

for i = step:dataCount-step
    n1 = dataSeries(i-step+1:i);
    n2 = dataSeries(i+1:i+step);
    mean1 = mean(n1);
    mean2 = mean(n2);
    c = (len1+len2)/(len1*len2);
    var1 = 1/len1*sum((n1 - mean1).^2);
    var2 = 1/len2*sum((n2 - mean2).^2);
    delta1 = len1*var1 + len2*var2;
    delta = delta1/(len1 + len2 - 2);
    t(i-step+1) = (mean1 - mean2)/sqrt(delta*c);
end
t
%%制图
figure(1)
plot(x,t,'r-','linewidth',1.5);
xlabel('t(year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimeNewRoman','FontSize',12);
axis([min(x),max(x),-4,4]);%-4,4,y轴坐标范围
hold on
plot(x,0*ones(i-step+1,1),'-.','linewidth',1);
plot(x,ttest*ones(i-step+1,1),':','linewidth',1);%更改线宽
plot(x,-ttest*ones(i-step+1,1),':','linewidth',1);
legend('t统计量','0.01显著水平');%这里的显著水平与前面的自由度v有关

1870-2019全球海表平均温度突变检验

clear all;close all;clc
%%     读取nc文件

lon = ncread('HadISST_sst.nc','longitude');
lan = ncread('HadISST_sst.nc','latitude');
sst = ncread('HadISST_sst.nc','sst');
time = ncread('HadISST_sst.nc','time');
time=time(1:1800);
sst(sst == -1000) = NaN;sst(sst < -1e30)=NaN;
%%      
time=time+datenum(1870,1,1);
time = double(time);
tstr=datestr(time,30);
mm=str2num(tstr(:,5:6));
msst=zeros(1800,1);
%求全球平均海表温度
for i=1:1800
    sst1 = sst(:,:,i);
    msst(i,1) = nanmean(nanmean(sst1,1),2);
end
%%    画全球平均海表温度的时间序列
figure(1)
set(gcf,'color','w')
plot(time,msst)
xlabel('时间（年）');
ylabel('温度（\circC）');
title('1870.01-2019.12全球平均海表温度时间序列曲线');
datetick('x','yyyy')

%%      滑动t检验
step = 16;n1 = 16;n2 = 16;%设置子序列步长
v = step+step-2;%自由度
ttest=1.70;%t检验值  alpha=0.1
k=1;
for j=1:300:1501
    time1=time(j:j+299,1);
    msst1=msst(j:j+299,1);
    k=k+1;
    t=[];
    for i = 1:length(time1)-step-step+1
         x1 = msst1(i:i+step-1);
         x2 = msst1(i+step:i+step+step-1);
         meanx1 = mean(x1);
         meanx2 = mean(x2);
         a = meanx1-meanx2;
         b = (n1+n2)/(n1*n2);
         varx1 = var(x1);
         varx2 = var(x2);
         c= ((n1-1)*varx1+(n2-1)*varx2)/(n1+n2-2);
         t1=a/sqrt(c*b);
         t=cat(1,t,t1);
    end
         x = time1(step:length(time1)-step);
         time1=time1+datenum(1870,1,1);
         time1 = double(time1);
         tstr=datestr(time1,30);
         
         figure(k)
         set(gcf,'color','w')
         plot(x,t,'r')
         axis([min(x),max(x),min(t),max(t)]);
         xlabel('year');
         ylabel('t统计量');
         datetick('x','yyyy')
         hold on
         x1=x(:,1);y1=ttest*ones(size(x1));
         plot(x1,y1,'b--')
         x2=x(:,1);y2=-ttest*ones(size(x2));
         plot(x2,y2,'b--')
         legend('t统计量','0.1显著水平');
end