clear all;
close all;

load precip_sno.txt;

  % first three values are year,month, and day
  % datenum in this instance converts time into a single 
  % numerical value so that we know
  % that consecutive days are separated by a value of 1
  dates = datenum(precip_sno(:,1:3));
  data(:,1:7) = precip_sno(:,4:10);
  %data(:,1)=ben lomond peak
  %data(:,2)=ben lomond trail
  %data(:,3)=tony grove lake
  %data(:,4)=snowbird
  %data(:,5)=thaynes canyon
  %data(:,6)=timpanogos
  %data(:,7)=payson
  
  label = ['BLP';'BLT';'TGL';'SNW';'THC';'TIM';'PAY'];
  
  [nd,nn] = size(dates);

  %sum data over each season
  dateold = 0;
  total = zeros(18,7);
  iy=0;
  for i=1:nd
   %increment season by looking for break between may and oct
    if(dates(i)-dateold > 100)
              iy=iy+1;
              yr(iy,1)=iy+1989;
    end
     total(iy,:)=total(iy,:)+data(i,1:7);  
     dateold = dates(i);
  end
  
  % plot blt and blp totals on one graph by year
  figure(1)
  plot(yr,total(:,1),'r',yr,total(:,2),'b')
  labelp = [label(1,:);label(2,:)];
  h=legend(labelp);
 
  %for convenience and to match notes, define blt as x and blp as y
  
x = total(:,2);
y = total(:,1);
n = 18;

xm = mean(x);
ym = mean(y);
xs = std(x,1);
ys = std(y,1);

% anomalies
xprime = x - xm;
yprime = y - ym;

% standardized anomalies

xstar = xprime/xs;
ystar = yprime/ys;

 % plot blt and blp standardized anomalies on one graph by year
  figure(2)
  plot(yr,ystar,'r',yr,xstar,'b')
  h=legend(labelp);

% brute force way to compute regression
% method 1

covar_1 = 0;
varx_1 = 0;
vary_1 = 0;

for i=1:n
covar_1 = covar_1 + xprime(i)*yprime(i);
varx_1 = varx_1 + xprime(i)*xprime(i);
vary_1 = vary_1 + yprime(i)*yprime(i);
end

covar_1 = covar_1/n;
varx_1 = varx_1/n;
vary_1 = vary_1/n;

b_1 = covar_1/varx_1;
r_1 = covar_1/sqrt(varx_1*vary_1);
sdx_1 = sqrt(varx_1);
sdy_1 = sqrt(vary_1);

% linear algebra way to compute

covar_2 = transpose(xprime)*yprime;
varx_2 = transpose(xprime)*xprime;
vary_2 = transpose(yprime)*yprime;
covar_2 = covar_2/n;
varx_2 = varx_2/n;
vary_2 = vary_2/n;

b_2 = covar_2/varx_2;
r_2 = covar_2/sqrt(varx_2*vary_2);
sdx_2 = sqrt(varx_2);
sdy_2 = sqrt(vary_2);

% compute estimate of y over range of x's 
%(we can choose whatever range of x's we want)
xhat = [ -15:1:15];
yhat = b_2 * xhat;


% add mean back in 
XH = xm + xhat;
YH = ym + yhat;

% plot estimate of y as line and plot original data as points
figure(3)
plot(XH,YH,'r-',x,y,'+','MarkerSize',10);
xlabel('Ben Lomond Trail','FontSize',16)
ylabel('Ben Lomond Peak','FontSize',16)

% scatter plot in terms of standardized anomalies
figure(4)
plot(xhat/xs,yhat/ys,'-',xstar,ystar,'+','MarkerSize',10);
xlabel('Ben Lomond Trail','FontSize',16)
ylabel('Ben Lomond Peak','FontSize',16)



% demo polytool
polytool(x,y)

% use matlab correlation function
rhop = corr(x,y,'type','Pearson');
rhos = corr(x,y,'type','Spearman');


%now start over for january only and estimate tgl from pay 
  
  dateold = 0;
  totalj = zeros(18,7);
  ij=0;
  for i=1:nd
    if(precip_sno(i,2)== 1)
     ij = precip_sno(i,1) - 1989;
      totalj(ij,:)=totalj(ij,:)+data(i,1:7);  
    end
  end
 
   % plot tgl and pay january totals on one graph by year
  figure(5)
  plot(yr,totalj(:,3),'r',yr,totalj(:,7),'b')
   labelp = [label(3,:);label(7,:)];
  h=legend(labelp);
  
  % going to try to estimate tgl from pay
  xj = totalj(:,7);
  yj = totalj(:,3);

% now another way of doing it in matlab is using the mldivide function 
% matrix left divide which solves a set of linear equations

xjm = mean(xj);
yjm = mean(yj);
xjs = std(xj,1);
yjs = std(yj,1);

xjp = xj - xjm;
yjp = yj - yjm;

xjstar = xjp/xjs;
yjstar = yjp/yjs;

% plot pay and tgl standardized anomalies on one graph by year
  figure(6)
  plot(yr,yjstar,'r',yr,xjstar,'b')
  h=legend(labelp);

%define a smaller range since not as much precip in january alone
xjhat = [ -2.5:.5:4];

%compute regression coefficient and correlation coefficient
bj = xjp\yjp;
rj = bj*xjs/yjs;

%estimate values over range
yjhat = xjhat*bj;


% add mean back in 
XJH = xjm + xjhat;
YJH = yjm + yjhat;

% plot january and linear fit 
figure(7)
plot(XJH,YJH,'r-',xj,yj,'+','MarkerSize',10);
xlabel('Payson','FontSize',16)
ylabel('Tony Grove','FontSize',16)



% figure in terms of standardized anomalies
figure(8)
plot(xjhat/xjs,yjhat/xjs,'-',xjp/xjs,yjp/xjs,'+','MarkerSize',10);
xlabel('Payson','FontSize',16)
ylabel('Tony Grove','FontSize',16)


% use matlab correlation function
rhojp = corr(xj,yj,'type','Pearson');
rhojs = corr(xj,yj,'type','Spearman');


