hw3.m

% BEGIN of HW3
% psilotum 2010/03/13
% for GNU Octave 

%% 1. Compute the mean and SE(mean) for the fish and copepod densities respectively,
% using both normal theory and non-parametric bootstrap

% import data first, and select the fish density
% and copepod density
enviANDdensity=xlsread("enviANDdensity.csv");
fish_and_copepod=enviANDdensity(:,11:12);
dfish=fish_and_copepod(:,1);
dcopepod=fish_and_copepod(:,2);

% (a) Normal theory
% calculate the mean 
mean_fc=mean(fish_and_copepod)

% calculate the standard error 
n=length(fish_and_copepod);
stderr_fc=sqrt(sum((fish_and_copepod-repmat(mean_fc,n,1)).^2)/n*(1/(n-1)))

% (b) Non-parametric boostrap

% call mybootstrap(obj,n,i) function, obj is the object
% n is the length of sample, and i is number of repeats
% repeat 1000 bootstrap
i=1000;
% calculate average of 1000 bootstrap means
dfish_bootstrap_means=mean(mybootstrap(dfish,n,i));
dcopepod_bootstrap_means=mean(mybootstrap(dcopepod,n,i));

dfish_bootstrap_avg_mean=mean(dfish_bootstrap_means)
dcopepod_bootstrap_avg_mean=mean(dcopepod_bootstrap_means)

% export results to file
save dfish_bootstrap_means.txt dfish_bootstrap_means;
save dcopepod_bootstrap_means.txt dcopepod_bootstrap_means;

% draw bootstrapped means with bootstrap 1000 times (with 100 bins)
% and export to pdf file
hist(dfish_bootstrap_means,100)
xlabel('Bootstrapped means of fish')
ylabel('Frequency')
print('figure1a.pdf', '-dpdf')
hist(dcopepod_bootstrap_means,100)
xlabel('Bootstrapped means of copepod')
ylabel('Frequency')
print('figure1b.pdf', '-dpdf')


% calculate the standard errors of bootstrap
stderr_dfish_bootstrap_means=sqrt(1/(i-1)* \
    sum((dfish_bootstrap_means-repmat(dfish_bootstrap_avg_mean,1,i)).^2))
stderr_dcopepod_bootstrap_means=sqrt(1/(i-1)* \
    sum((dcopepod_bootstrap_means-repmat(dcopepod_bootstrap_avg_mean,1,i)).^2))

%-----------------------------------------------------
%% 2. Compute the median and bootstrapped SE(median) for the fish 
%  and copepod densities

median_fc=median(fish_and_copepod)

% calculate the median of bootstrapped densities
bstrp_f=mybootstrap(dfish,n,i);
bstrp_c=mybootstrap(dcopepod,n,i);

% median of bootstrapped fish & copepod densities
fxi=median(bstrp_f);
cxi=median(bstrp_c); 
fxbar=mean(fxi) % average of bootstrapped median of fish
cxbar=mean(cxi); % average of bootstrapped median of copepod

% export results to files
save median_bootstrap_fish.txt fxi;
save median_bootstrap_copepod.txt cxi;

% calculate the standard error of bootstrapped median
se_dfish_bootstrap_median=sqrt(1/(i-1)*sum((fxi-fxbar).^2))
se_copepod_bootstrap_median=sqrt(1/(i-1)*sum((cxi-cxbar).^2))

% plot the histograms  with bootstrap 1000 times
hist(fxi,30)
xlabel('Bootstrapped medians of fish')
ylabel('Frequency')
print('figure2a.pdf', '-dpdf')
hist(cxi,30)
xlabel('Bootstrapped medians of copepod')
ylabel('Frequency')
print('figure2b.pdf', '-dpdf')

%-----------------------------------------------------
%% 3.1 Plot fish (Y) v.s copepod (X) and the regression 
%  line.

% find beta_0 and beta_1
solve=[ones(length(dcopepod), 1) dcopepod]\dfish
beta0=solve(1,:);
beta1=solve(2,:);

% plot points
plot(dcopepod,dfish,'*')
hold on
    yfit=beta0+beta1*dcopepod;
    plot(dcopepod,yfit,'r')
hold off

xlabel('Copepod densities (ind./m^3)')
ylabel('Fish densities (ind.)')
print('figure3a.pdf', '-dpdf')


%% 3.2 Compute the regression coefficients for
%  fish = beta0+beta1*copepod and bootstrapped SE(beta0) and SE(beta1). Plot 
%  the histogram of bootstrapped beta0 and beta1 with bootstrap 1000 times.

% find the beta_0 and beta_1
bstrp_fc=mybootstrap(fish_and_copepod,n,i); 

bootstrap_solve=[];
for k=1:i
    bootstrap_solve=[bootstrap_solve, [ones(length(dcopepod),1) bstrp_fc(:,k,2)]\bstrp_fc(:,k,1)];
end
 
bootstrap_b0=bootstrap_solve(1,:);
bootstrap_b1=bootstrap_solve(2,:);
 
save bootstrap_b0.txt bootstrap_b0
save bootstrap_b1.txt bootstrap_b1
  
% calculate the standard error of beta0 and beta1
se_beta0=sqrt(1/(i-1)* \
     sum((bootstrap_b0-repmat(mean(bootstrap_b0),1,i)).^2))
  
se_beta1=sqrt(1/(i-1)* \
     sum((bootstrap_b1-repmat(mean(bootstrap_b1),1,i)).^2))
  
hist(bootstrap_b0,100)
xlabel('Bootstrapped coefficient beta_0')
ylabel('Frequency')
print('figure3b.pdf', '-dpdf')
hist(bootstrap_b1,100)
xlabel('Bootstrapped coefficient beta_1')
ylabel('Frequency')
print('figure3c.pdf', '-dpdf')