%  These are a few commands that may help with homework 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Problem 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Load the appropriate data
load madison_precip.txt

%  The first column is the year, the second is month, the third is day,
%  and the last is measured precip.  I have replaced all "trace" days
%  with 0.005.  Let's extract the data:
yr = madison_precip(:,1);
mon = madison_precip(:,2);
day = madison_precip(:,3);
prec = madison_precip(:,4);

%  Now, let's just get May Precipitation
may_indices = find(mon == 5);
pmay = prec(may_indices);

%  We can see how many days we have in our record by:
ndaytot = length(pmay)

%  Let's suppose we'd like to find how many days have a measured 0.01
%  inch of precipitation in May:  (leave the semicolon off the end, so we
%  can see the output)
nday0p01 = sum(pmay == 0.01)

figure(1); orient tall; wysiwyg
subplot(2,1,1);
hist(nday0p01, 0:31);

%  Note that we could use >, <, >=, <=, or the compliment ~ (not) by ~=
%  (not equal)

%  Suppose we'd like to identify how many instances we have at least two
%  consecutive days of rain IN MAY:

may_indices1 = find((mon == 5) & (day <= 30));  %  All days but the last
                                                %  day in May
may_indices2 = may_indices1 + 1;                %  The next day

%  Now, find consecutive days in May by finding when both day 1 AND day 2
%  have precipitation:
ntwodaysp = sum((prec(may_indices1) > 0) & (prec(may_indices2) > 0))



%  Now, let's look at the number of days of rain in a given month.  Start
%  by reformatting the data:
pmay = reshape(pmay, 31, 48);  %  Reshape so that each row represents a
                               %  different day in May, and every column
                               %  is a different year
%  Note that the above command is the same as:
%  for i = 1:48;  %  there are 48 years
%    tem(:,i) = pmay(31*(i-1)+[1:31]);
%  end
%  pmay = tem;

raindays = sum(pmay > 0);



%  Now, you can calculate the theoretical binomial distribution via the
%  'binopdf' function (type 'help binopdf'):
p = 0.3642;  %  'p' is the probability that there will be rain on any
             %  given day.  How did I calculate that?
fx = binopdf(0:31, 31, p);

  
%%% NOTE:  You will need the statistics toolbox to use 'binopdf' or
%%% 'binocdf'.  

  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Problem 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  To load madison temperature, do the following:
load madison_temperature.txt
mt = madison_temperature;
yr = mt(:,1);
mon = mt(:,2);
day = mt(:,3);
meant = mt(:,4);

%  Now, get January only temperature
jant = meant(mon == 1);

%  Suppose we'd like to identify the quartiles for jant:
[jant_sort, jant_ind] = sort(jant);
q1 = jant_sort(floor(length(jant)/4));
med = median(jant_sort);
q3 = jant_sort(ceil(3*length(jant)/4));




%  Make a probability plot.  

%  For this exercise, we'll make a probability plot for a normally
%  distributed random variablle.  First, let's define the random
%  variable: 
y = 6*randn(100,1) + 30;

%  You should try this with other random distributions:
%  y = gamrnd(5, 2, 1000, 1);
%  y = binornd(30, 0.3, 1000, 1);

%  To make a probability plot, we'd like to scale the x-axis so that a
%  normal distribution is a straight line.  Let's define the x-axis.
%  What we'd like to do is divide the normal distribution into segments
%  of equal area, with each 'x' value centered on each bin.  We can
%  achieve this by noting that equal area implies equal probability, so
%  we just take the inverse of the cumulative normal distribution:
%
prob = ([1:length(y)]-0.5)/length(y);  %  The center of each bin
x = norminv(prob, 0, 1);  %  Plot this against a standard normal
                          %  distribution

%  Now, get a theoretical normal distribution with the same parameters as
%  'y'.  
yt = norminv(prob, mean(y), std(y));

%  Now, plot the results:
figure(1); clf
orient tall; wysiwyg;

subplot(2,1,1);

%  Plot the empirical and theoretical distributions on a graph:
  pp = plot(x, sort(y), '.k', x, yt, '-k');

%  And make it look pretty
  axis square
  xlim([-3.5 3.5]);

%  Now, let's label the x-axis so that it represents the cumulative
%  probability:
  xtic_prob = [.001 0.01 0.05 .1 .2 .5 .8 .9 .95 .99 .999];
  set(gca, 'XTick', norminv(xtic_prob), 'XTickLabel', xtic_prob);
  
%  Make it pretty, again:
  grid on
  xl = xlabel('Cumulative Probability', 'fontsize', 12);
  yl = ylabel('NORM Random Variable', 'fontsize', 12);
  t1 = title('Probability Plot for a Normal Random Variable', 'fontsize', 14);