%% MATLAB Recipes for Earth Sciences - Chapter 3

%% Section 3.3
clear

%%
corg = load('organicmatter_one.txt');

%%
plot(corg,zeros(1,length(corg)),'o')

%%
hist(corg)

%%
[n,v] = hist(corg)

%%
hist(corg,8)

%%
min(corg) + 0.5*range(corg)/8
max(corg) - 0.5*range(corg)/8

%%
range(corg)/8

%%
v = 9.75 : 0.65 : 14.30;

%%
hist(corg,v)

%%
n = hist(corg,v)

%%
mean(corg)

%% 
median(corg)

%%
quantile(corg,[.25 .50 .75])

%%
v(find(n == max(n)))

%%
v(n == max(n))

%%
boxplot(corg)

%%
var(corg)

%%
std(corg)

%%
skewness(corg)

%%
kurtosis(corg)

%%
corg(25,1) = NaN;

%%
mean(corg)

%%
nanmean(corg)

%%
clear

%%
sodium = load('sodiumcontent_one.txt');

%%
hist(sodium,11)

%%
[n,v] = hist(sodium,11)

%%
mean(sodium)

%%
median(sodium)

%%
v(find(n == max(n)))

%%
skewness(sodium)

%%
clear

%%
sodium = load('sodiumcontent_two.txt');

%%
[n,v] = hist(sodium,7);

%%
mean(sodium)

%%
median(sodium)

%%
v(find(n == max(n)))

%%
sodium(51,1) = 1.5;

%%
hist(sodium,7)

%%
mean(sodium)

%%
median(sodium)

%%
v(find(n == max(n)))

%% Section 3.5
clear

%%
x = 9 : 0.1 : 15;
pdf = normpdf(x,12.3448,1.1660);
cdf = normcdf(x,12.3448,1.1660);
plot(x,pdf,x,cdf)

%% Section 3.6
clear

%%
load('organicmatter_two.mat');

%%
[n1,x1] = hist(corg1);
[n2,x2] = hist(corg2);

%%
h1 = bar(x1,n1);
hold on
h2 = bar(x2,n2);
hold off

set(h1,'FaceColor','none','EdgeColor','r')
set(h2,'FaceColor','none','EdgeColor','b')

%%
[n1,x1] = hist(corg1);
[n2,x2] = hist(corg2);

%%
na = length(corg1);
nb = length(corg2);
ma = mean(corg1);
mb = mean(corg2);
sa = std(corg1);
sb = std(corg2);

%%
tcalc = abs((ma-mb))/sqrt(((na+nb)/(na*nb)) * ...
    (((na-1)*sa^2+(nb-1)*sb^2)/(na+nb-2)))

%%
tcrit = tinv(1-0.05/2,na+nb-2)

%%
[h,p,ci,stats] = ttest2(corg1,corg2,0.05)

%%
clear

%%
load('organicmatter_three.mat');

%%
[n1,x1] = hist(corg1);
[n2,x2] = hist(corg2);

%%
na = length(corg1);
nb = length(corg2);
ma = mean(corg1);
mb = mean(corg2);
sa = std(corg1);
sb = std(corg2);

%%
tcalc = abs((ma-mb))/sqrt(((na+nb)/(na*nb)) * ...
    (((na-1)*sa^2+(nb-1)*sb^2)/(na+nb-2)))

%%
tcrit = tinv(1-0.05/2,na+nb-2)

%%
[h,p,ci,stats] = ttest2(corg1,corg2,0.05)

%% Section 3.7
clear

%%
load('organicmatter_four.mat');

%%
s1 = std(corg1)
s2 = std(corg2)

%%
df1 = length(corg1) - 1
df2 = length(corg2) - 1

%%
if s1 > s2
slarger = s1
ssmaller = s2
else
slarger = s2
ssmaller = s1
end

%%
Fcalc = slarger^2 / ssmaller^2
Fcrit = finv(1-0.05/2,df1,df2)

%%
[h,p,ci,stats] = vartest2(corg1,corg2,0.05)

%%
clear

%%
load('organicmatter_five.mat');

%%
s1 = std(corg1);
s2 = std(corg2);

%%
df1 = length(corg1) - 1;
df2 = length(corg2) - 1;

%%
if s1 > s2
slarger = s1;
ssmaller = s2;
else
slarger = s2;
ssmaller = s1;
end

%%
Fcalc = slarger^2 / ssmaller^2
Fcrit = finv(1-0.05/2,df1,df2)

%%
[h,p,ci,stats] = vartest2(corg1,corg2,0.05)

%% Section 3.8
clear

%%
corg = load('organicmatter_one.txt');

%%
v = 9.40 : 0.74 : 14.58;
n_obs = hist(corg,v);

%%
n_exp = normpdf(v,mean(corg),std(corg));

%%
n_exp = n_exp / sum(n_exp);
n_exp = sum(n_obs) * n_exp;

%%
subplot(1,2,1), bar(v,n_obs,'r')
subplot(1,2,2), bar(v,n_exp,'b')

%%
chi2calc = sum((n_obs - n_exp).^2 ./ n_exp)

%%
chi2crit = chi2inv(1-0.05,5)

%%
[h,p,stats] = chi2gof(corg)

%% Section 3.9
clear

%%
randn('seed',0)
data = 6.4 + 1.4*randn(100,1);

%%
v = 2 : 10;
hist(data,v)
n = hist(data,v);

%%
[muhat,sigmahat] = normfit(data)

%%
x = 2 : 1/20 : 10;
y = normpdf(x,muhat,sigmahat);
y = trapz(v,n) * y/trapz(x,y);
bar(v,n), hold on, plot(x,y,'r'), hold off

%%
phat = mle(data,'distribution','normal');

%%
x = 2 : 1/20 : 10;
y = normpdf(x,phat(:,1),phat(:,2));
y = trapz(v,n) * y/trapz(x,y);

%%
bar(v,n), hold on, plot(x,y,'r'), hold off

%%
clear
randn('seed',0)
ya = 6.4 + 1.4*randn(100,1);
yb = 13.3 + 1.8*randn(100,1);
data = vertcat(ya,yb);

%%
v = 0 : 30;
hist(data,v)
n = hist(data,v);

%%
gmfit = gmdistribution.fit(data,2)

%%
x = 0 : 1/30 : 20;
y1 = normpdf(x,gmfit.mu(1,1),gmfit.Sigma(:,:,1));
y2 = normpdf(x,gmfit.mu(2,1),gmfit.Sigma(:,:,2));

%%
y1 = gmfit.PComponents(1,1) * y1/trapz(x,y1);
y2 = gmfit.PComponents(1,2) * y2/trapz(x,y2);

%%
y = y1 + y2;
y = trapz(v,n) * y/trapz(x(1:10:end),y(1:10:end));

%%
bar(v,n), hold on, plot(x,y,'r'), hold off


