%% Exercise 2: Trends in ODP 659 Benthic Oxygen Isotope Record
% Loading data
clear
data = load('odp659data.txt');

%%
% Plotting again like yesterday
figure1 = figure('Color',[1 1 1],...
    'Position',[150 400 550 260]);
axes1 = axes('Box','on',...
    'XDir','reverse','YDir','reverse',...
    'Units','Centimeters','Position',[2 2 16 6],...
    'LineWidth',0.6,'FontName','Helvetica','FontSize',12);
line1 = line(data(:,1),data(:,2));
set(line1,'Color',[0.8 0.3 0.1])
set(get(axes1,'Title'),...
    'String','ODP 659 Benthic Oxygen Isotope Record',...
    'FontName','Helvetica','FontSize',12)
set(get(axes1,'XLabel'),...
    'String','Age- Thousands of Years Before the Present',...
    'FontName','Helvetica','FontSize',12)
set(get(axes1,'YLabel'),...
    'String','d18O- Permille v. PDB',...
    'FontName','Helvetica','FontSize',12)

%%
% Three sections separated at 1,000 kyrs and 2,700 kyrs BP
part1 = data(1:251,:);
part2 = data(252:669,:);
part3 = data(670:1194,:);

%%
% Plotting the three parts
subplot(3,1,1), line(part1(:,1),part1(:,2),'Color',[0.8 0.3 0.1])
subplot(3,1,2), line(part2(:,1),part2(:,2),'Color',[0.8 0.3 0.1])
subplot(3,1,3), line(part3(:,1),part3(:,2),'Color',[0.8 0.3 0.1])

%%
% Histogram[h,p,ci,stats] = ttest2(part1(:,2),part2(:,2),0.05)s of the three parts, seems to be Gaussian
subplot(1,3,1), hist(part1(:,2))
subplot(1,3,2), hist(part2(:,2))
subplot(1,3,3), hist(part3(:,2))

%%
% Statistics of the three parts assuming that the data are Gaussian
statistics(1,1) = mean(part1(:,2)); statistics(1,2) = var(part1(:,2));
statistics(2,1) = mean(part2(:,2)); statistics(2,2) = var(part2(:,2));
statistics(3,1) = mean(part3(:,2)); statistics(3,2) = var(part3(:,2));
statistics

%%
% T-test assuming that the data are Gaussian. H=1 indicates that the null
% hypothesis can be rejected at the 5% level in both cases, so the diff in
% the means is significant. Same with F-test.
[h,p,ci,stats] = ttest2(part1(:,2),part2(:,2),0.05)
[h,p,ci,stats] = ttest2(part2(:,2),part3(:,2),0.05)
[h,p,ci,stats] = vartest2(part1(:,2),part2(:,2),0.05)
[h,p,ci,stats] = vartest2(part2(:,2),part3(:,2),0.05)

%%
% Polynomial fits
p11 = polyfit(part1(:,1),part1(:,2),1);
p22 = polyfit(part2(:,1),part2(:,2),1);
p33 = polyfit(part3(:,1),part3(:,2),1);

%%
[r1,p1] = corrcoef(part1(:,1),part1(:,2))
tcalc = r1(2,1) * ((length(part1(:,1))-2)/(1-r1(2,1)^2))^0.5
tcrit = tinv(0.95,length(part1(:,1))-2)
[r2,p2] = corrcoef(part2(:,1),part2(:,2))
tcalc = r2(2,1) * ((length(part1(:,1))-2)/(1-r1(2,1)^2))^0.5
tcrit = tinv(0.95,length(part1(:,1))-2)
[r3,p3] = corrcoef(part3(:,1),part3(:,2))
tcalc = r3(2,1) * ((length(part1(:,1))-2)/(1-r1(2,1)^2))^0.5
tcrit = tinv(0.95,length(part1(:,1))-2)

%%
% Plotting
subplot(3,1,1) 
plot(part1(:,1),part1(:,2),'o'), hold on
plot(part1(:,1),polyval(p11,part1(:,1)),'r'), hold off
subplot(3,1,2) 
plot(part2(:,1),part2(:,2),'o'), hold on
plot(part2(:,1),polyval(p22,part2(:,1)),'r'), hold off
subplot(3,1,3) 
plot(part3(:,1),part3(:,2),'o'), hold on
plot(part3(:,1),polyval(p33,part3(:,1)),'r'), hold off

%%
% Alternatives? How do find the breakpoints in the records?
% (1) Mann-Whitney-Test of sliding windows (by M Trauth), used in Trauth et al. QSR 2009
% (2) Breakfit Regression (by M Mudelsee), used in Trauth et al. QSR 2009
% (3) Bayesian Regression (by N Schuetz), not yet published but in press Physical Review E
% (4) Recurrence Networks (by J Donges), accepted for pub in Nonlin Proc Geosphysics
% (5) Rampfit Method (by M Mudelsee)