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


close all;
clear all;
clc;
clf;

%% In this new version we consider three n and try to find the correct scaling. 
m=1000;
n1=10;
n2=100;
n3=1000;
nMax=max([n1,n2,n3]);
%% Generate data.

%%(1)Uniform distribution
%Data=rand(m,nMax);

%%(2) Exponential distribution
%Lambda=1;
%Data=exprnd(Lambda,m,nMax);

%%(3) Pareto distribution. Alpha and Xm are paramters in your lecture notes.
%%K,Sigma and Theta are parameters used by matlab. We do a little
%%translation first.
  Alpha=4;Xm=1;
  K=1/Alpha;Sigma=Xm/Alpha;Theta=Xm;
  Data=gprnd(K,Sigma,Theta,[m,nMax]);

%% Take the max of each row.
maxData1=zeros(m,1);
maxData2=zeros(m,1);
maxData3=zeros(m,1);
for i=1:m
    maxData1(i)=max(Data(i,1:n1));
    maxData2(i)=max(Data(i,1:n2));
    maxData3(i)=max(Data(i,1:n3));
end

%% To see data collapse we need to scale extreme values M by (M-b)/a

%%(1)Uniform
%a=[1/n1,1/n2,1/n3];
%b=[1,1,1];

%%(2)Exponential
%a=[1,1,1];
%b=[log(n1)/Lambda,log(n2)/Lambda,log(n3)/Lambda];

%%(3)Pareto
a=[Xm*(n1^K),Xm*(n2^K),Xm*(n3^K)];
b=[0,0,0];





%% Plot
hold on
 [f1,x1]=ksdensity((maxData1-b(1))./a(1)); %Estimate the pdf of scaled data
 plot(x1,f1,'og','Markersize',15);
 
 [f2,x2]=ksdensity((maxData2-b(2))./a(2));
 plot(x2,f2,'xr','Markersize',15);

 [f3,x3]=ksdensity((maxData3-b(3))./a(3));
 plot(x3,f3,'sb','Markersize',15);

 %%Theoratical curve. Here we fit parameters with gevfit. You can also try
 %%to compute the actual parameters and plot a better theory curve.
 param=gevfit((maxData3-b(3))./a(3)); 
 xTheory=x3;
 yTheory=gevpdf(xTheory,param(1),param(2),param(3));
 plot(xTheory,yTheory,'k');
 
xlabel('x');
ylabel('PDF');
legend('n=10','n=100','n=1000','Theory');
hold off