%% Class for CO907, Fri, Week 1.
%__________________________________________
%Author: Jiarui Cao
%Date: 8 Nov 2013
%Purpose: Sheet1. Queation 1.7a


close all;
clear all;
clc;
clf;

rng('shuffle') ;
m=500;

%% Craet a m*1000 random matrix as our data.

%Data=exprnd(1,m,1000); % exponential distribution. 
%Data=rand([m,1000]); % uniform distribution.
Data=gprnd(0.1,1,1,[m,1000]); %generalized pareto distribution.

%% Compute the mean and deviation of ALL data.
DataMean=mean(mean(Data));
DataStd = std2(Data); %std2 returns std deviation of a matrix


%% Creat 3 sum arrary
s10=zeros(m,1);
s100=zeros(m,1);
s1000=zeros(m,1);
n1=10;
n2=100;
n3=1000;

%% Main loop. For each i, we deal with the ith row of our data matrix.
for i=1:m
    s10(i)=sum(Data(i,1:n1)); %sum the first n1 data.
    s10(i)= (s10(i)-(n1*DataMean))/(DataStd*sqrt(n1)); %scale it.
    
    s100(i)=sum(Data(i,1:n2));%sum the first n2 data.
    s100(i)= (s100(i)-(n2*DataMean))/(DataStd*sqrt(n2));
    
    s1000(i)=sum(Data(i,1:n3));%sum the first n3 data.
    s1000(i)= (s1000(i)-(n3*DataMean))/(DataStd*sqrt(n3));
end


%% Plot the probability density estimate.
hold on;
[f1,x1]=ksdensity(s10);% look up 'ksdensity' in matlab help file if necessary
plot(x1,f1,'or'); % 'or' means use symbol 'circle' and color 'red';

[f2,x2]=ksdensity(s100);
plot(x2,f2,'xb'); % cross and blue;

[f3,x3]=ksdensity(s1000);
plot(x3,f3,'sg'); %square and green;


%% Plot pdf of a std normal distribution  
xx=-4:0.01:4;
yy=normpdf(xx,0,1); 
plot(xx,yy,'k','linewidth',2);
hold off

title('PDF plots for CLT');
xlabel('x');
ylabel('PDF');
legend('S10','S100','S1000','Std Normal','Location','South');
set(gca,'Yscale','log'); %set gca (current axeis handle) to be y-log. Look up 'set' for more interesting functions.