MATLAB滑动T检验
⽤MATLAB做滑动T检验
滑动t检验是通过考察两组样本平均值的差异是否显著来检验突变。
基本思想是：把⼀⽓候序列中两段⼦序列均值有⽆显著差异看作来⾃两个总体均值有⽆显著差异的问题来检验。如果两段⼦序列的均值差异超过了⼀定的显著性⽔平，则可以认为有突变发⽣。
本篇博客中的程序1⽐较结构性⽐较差，⽐较乱，程序2的可读性更好
嘿嘿，第⼀个程序是我⾃⼰编的，有很⼤改进空间，第⼆个程序是⽼师给的，⽅便改参数，直接copy就能⽤，真⾹。
程序1
clear all;close all;clc
%% data read
data = xlsread('Tair.xlsx');
time  = data(:,1);
m_sst = data(:,2);
%%
step_1 =10;step_2 =10;%设置两组⼦序列的步长
length_1 = length(time);%变量的命名不能跟函数名相同，不然再次使⽤该函数时会报错
%⽤三个循环，⽅便每⼀个变量的查看，也可以将三个循环整合成⼀个
x1 =[];ss_1 =[];
%求第⼀个⼦序列的平均和⽅差并放在数组X1和SS_1中
for i =1:length_1 - step_1 - step_2
x1_bar = mean(m_sst(i:i+step_1-1));
a = m_sst(i:i+step_1-1);
s_1 =0;
for j =1:step_1
s_1 =  s_1 +(a(j)- x1_bar)^2;%可以直接⽤sum函数求和
end
s_1 = s_1./step_1;
x1 = cat(1,x1,x1_bar);
ss_1 = cat(1,ss_1,s_1);%使⽤cat函数时，不能将拼接的两个变量打乱，否则会得到⼀个顺序相反或者错乱的数组endfontweight属性bold
%求第⼆个⼦序列的平均和⽅差并放在数组X2和SS_2中
x2 =[];ss_2 =[];
for k =1+step_1:length_1 - step_2
x2_bar = mean(m_sst(k:k + step_2 -1));
a = m_sst(k:k + step_2 -1);
x2 = cat(1,x2,x2_bar);
s_2 =0;
for l =1:step_2
s_2 = s_2 +(a(l)- x2_bar)^2;
end
s_2 = s_2./step_2;
ss_2 = cat(1,ss_2,s_2);
end
%求滑动t值并拼接成⼀个矩阵
sss =[];
length_2 = length(x2);
tt =[];
for m =1:length_2
b = sqrt(((step_1 -1)*ss_1(m)+(step_2 -1)*ss_2(m))/(step_1+step_2 -2));
sss = cat(1,b,sss);
t =(x1(m)-x2(m))/(b*sqrt(1/step_1 +1/step_2));
tt = cat(1,tt,t);
end
length_3 = length(tt);
set(0,'defaultfigurecolor','w');%画布背景设置
set(gcf,'units','normalized','position',[0.250.250.50.5]);
plot(time(step_1 +1:length_3 + step_1),tt,'linewidth',2);
hold on
ttest =2.878;% significant level  alpha=0.01;
plot(time(step_1 +1:length_3 + step_1),zeros(length_3,1),'-.','linewidth',1);
plot(time(step_1 +1:length_3 + step_1),ttest*ones(length_3,1),':','linewidth',3);% line of the significant level
plot(time(step_1 +1:length_3 + step_1),-ttest*ones(length_3,1),':','linewidth',3);
H=legend('t','0.01 significant level');% explain
title('滑动t-检验检测1911-1995年中国年平均⽓温等级序列的突变','fontweight','bold','fontsize',15);
xlabel('t(year)','FontName','TimesNewRoman','FontSize',12,'fontweight','bold');
ylabel('Moving T Test','FontName','TimeNewRoman','FontSize',12,'fontweight','bold');
box off;
结果1
公式
啊，对了，附上滑动T检验的公式，写着写着就忘了
两个⼦样本的样本长度分别为：n1，n2
均值分别为：
⽅差分别为：S1,S2
两个滑动⼦序列
⽰意图
程序2
clear all;close all;clc
%%打开数据
f ='';
time = ncread(f,'time');
lat  = ncread(f,'latitude');
lon  = ncread(f,'longitude');
sst  = ncread(f,'sst');
%%将异常值变为NAN
sst(sst ==-1000)=NaN;sst(sst ==-1.8)=NaN;sst(sst <-1e30)= NaN; %%求全球海表平均温度
for t =1:size(sst,3)
sst_d    = squeeze(sst(:,:,t));
m_sst(t)= nanmean(sst_d(:));
end
%m_sst = squeeze(nanmean(squeeze(nanmean(sst,1)),1));
%%滑动T检验
step =10;% length of subsequence
v = step+step-2;% degreee of freedom
ttest =2.878;% sinnificant level  alpha=0.01;
len1 = step;
len1 = step;
len2 = step;
length_1 = length(m_sst);
x = time(step:length_1 - step);
%将时间序列x变为正常的时间序列
x = x + datenum(1870,01,01);
for i = step:length_1 - step
n1 = m_sst(i-step+1:i);
n2 = m_sst(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-1)*var1 +(len2-1)*var2;
delta = delta1/(len1 + len2 -2);
t(i-step+1)=(mean1 - mean2)/sqrt(delta*c);
end
%%趋势图
figure(1)
set(0,'defaultfigurecolor','w');%画布背景设置
set(gcf,'units','normalized','position',[0.10.10.80.8]);
subplot(2,1,1);
time = double(time);
time = time + datenum(1870,01,01);
plot(time,m_sst,'k-');
datetick('x','mm/dd/yyyy','keepticks');%改时间坐标轴
xlabel('t(year)','FontName','TimesNewRoman','FontSize',12,'fontweight','bold');
ylabel('SST','FontName','TimeNewRoman','FontSize',12,'fontweight','bold');
title('全球海表平均温度时间序列图','fontweight','bold','fontsize',20);
box off;
%%滑动T检验图
subplot(2,1,2);
plot(x,t,'b-','linewidth',1);
datetick('x','mm/dd/yyyy','keepticks');%改时间坐标轴
xlabel('t(year)','FontName','TimesNewRoman','FontSize',12,'fontweight','bold');
ylabel('Moving T Test','FontName','TimeNewRoman','FontSize',12,'fontweight','bold');
axis([min(x),max(x),-6,6]);% y axis limitation
hold on
plot(x,0*ones(i-step+1,1),'-.','linewidth',1);
plot(x,ttest*ones(i-step+1,1),':','linewidth',3);% line of the significant level
plot(x,-ttest*ones(i-step+1,1),':','linewidth',3);
H=legend('t','0.01 significant level');% explain
title('⽤滑动T检验检测给出1870.01-2019.12年全球海表平均温度序列的突变','fontweight','bold','fontsize',20); box off;
结果2