%% Temple university CIS 2033
% A modern Introduction to Probability and Statistics
% Chapter 7: Expectation and variance
%
% Author: Djordje Gligorijevic
% email: gligorijevic@temple.edu

%% North east region of United States precipitation
% We will briefly play with data on precipitation in north east region
% (where we are now!). 
% First load the data

load precipitation_north_east.mat;

% As these are real life data there are supposed to be a lots of missing
% values! The missing data have value NaN in the dataset. 
% What we wish to do is to preprocess this data before we can use it.
% There are many ways to preprocess the dataset with missing value, but
% what we are going to do is simple. We will just replace NaN value with
% the zeros.
% In the class you will be required to replace the missing values with the
% average precipitation for each of 100 stations in the dataset!
y(isnan(y)) = 0;

% Now lets analyze mean value of precipitation for the north east region of
% united states:
meanPrecipitation = mean(mean(y));

% Now lets see its variance:
varPrecipitation = var(var(y));

% Great, now let us plot all 100 stations on a single plot over time given:
figure;
hold on;
for i = 1:size(y,1)
   plot(y(i,:)); 
end
title('Precipitation ammount for each station');
ylabel('Precipitation ammount');
xlabel('Time');

% It's a little crowdy isn't it? What could we possibly do to improve the
% plot?
% There is a way, that using mean and variance of the data allows us to
% make a nicer plot! 
% In each time step, we can measure mean precipitation over all 100
% stations, then measure variance as well, and we can plot both onto a
% single plot.
meanPrecipiation = NaN(size(y,1));
varPrecipitation = NaN(size(y,1));

for i = 1:size(y,2)
   meanPrecipiation(i) = mean(y(:,i));
   varPrecipitation(i) = var(y(:,i));
end

figure;
hold on;
plot(meanPrecipiation(:), 'r', 'lineWidth',1.5);
plot(meanPrecipiation(:) + 1.96*sqrt(varPrecipitation(:)), 'b--', 'lineWidth', 0.8);
plot(meanPrecipiation(:) - 1.96*sqrt(varPrecipitation(:)), 'b--', 'lineWidth', 0.8);
title('Precipitation ammount for each station');
ylabel('Precipitation ammount');
xlabel('Time');