%% Problem 1
% As the integral of a function equals to the area of the x-y plane under
% the curve, I use random numbers method to find this area.

fprintf('\nProblem 1 Results\n\n')

%A vector of random numbers in x-y plane:
N=1000;
X=pi/2.*rand(N,1);
Y=rand(N,1);

% Calculate the total area of the x-y plane in special region which is
% defined by lower and upper bound of integral
Area_Rectangular=1*pi/2;

% Define the function
f=@(x) cos(x);

% Plot the random numbers
plot(X,Y,'.')
xlabel('X')
ylabel('Y')
axis equal
hold on

% Plot the function
plot(X,f(X),'r.')
legend('Random Numbers','Cos(X)')
hold off

% find the numbers of dots under the curve
N_F=sum(Y<f(X));

% Calculate the area, I know that the rate of given area in order of the
% total area is obtained by N_F/N
% I expect this result should be close to 1 from the known solution
Area=N_F/N*Area_Rectangular;

% display the result
fprintf('The integral of cos(x) from 0 to pi/2 is %g\n', Area)

%% Problem 2

fprintf('\nProblem 2 Results\n\n')

k_perday=1;      % Rate constant in days
T_days=4;        % Total time in days. 

% I want to simulate this process with 'small' dt values.  'small'
% means when P=k*dt << 1.  For conceptual convienience, lets convert to
% minutes. 

MPD=24*60;  %Minutes per day. 

dt =1 ;             %min
k_min=k_perday/MPD; %min^-1
P=k_min*dt;         %Probabilty of an event in each interval

%P is <<1 so the size of dt looks appropriate. 

%Number of 1 min intervals in 4 days. 
N_intervals = 4*MPD;    

%Check that the expected counts over the total time period is correct:
lambda=k_min*N_intervals*dt

%Simulate experiment 1000 times. 
Repeats = 10000;        

%Simulate 1000 experiments
Events=double(rand(N_intervals,Repeats)<P);

%find number of events for each repeat trial
Counts=sum(Events);

%Show the distribution
figure
hist(Counts,(0:10))
xlabel('Time')
ylabel('Counts')
s=sprintf('Lambda = %g',lambda)
title(s)

%The probability that counts = 0 
P_0 = sum(Counts==0)/Repeats

fprintf('\n Estimated probability of no events in 4 days from simulations: %g\n\n',P_0)

%Check against the Poisson PMF
P_0PMF=poisspdf(0,4)

fprintf('\n Probability of no events in 4 days from PMF: %g\n\n',P_0PMF)


%% Problem 3
% In this problem, I like to compare the exact value of exponential
% function with the given limitation relation

fprintf('\nProblem 3 Results\n\n')

% Give value to x 
x=1;

% Find the exact value for x
Y_exact=exp(x);

% to take the n goes to inifinity, I need to define a 'for' loop and give n
% as a big value.
N=1000000;

% find the relation to find the exponential value
for i=1:N
    n(i)=i;
    Y(i)=(1+x/n(i))^n(i);
    div(i)=abs(Y(n(i))-Y_exact);
end

% Plot the result
figure
loglog(n,div)
title('Deviation of exp(1) from (1+1/n)^n')
xlabel('log(n)')
ylabel('Error')
legend('Deviation')
%% Problem 4
% Look at the given link in the problem. I need to create poisson
% distribution for special lambdas.
% Here is two methods to solve this problem,

fprintf('\nProblem 4 Results\n\n')

% First method
% Give some values for lambda
lambda=[1 4 10];


% Here is a solution using poissrnd: 

%Repeat draws to estimate the distribution
N=100000;   
Kvec=(0:20);
figure
hold on
LineColor(1,:)=[1 0.5 .2];
LineColor(2,:)=[0.5 0 .5];
LineColor(3,:)=[0 0 1];

for ll=1:3 %Loop over lambda
    L=lambda(ll);
   
    %Draw random numbers from Poisson:
    K=poissrnd(ones(N,1)*L);
    
    %Use histogram to count how many are in each integer
    [Counts]=hist(K,Kvec);
    
    %Normalize:
    Counts_norm=Counts/sum(Counts);
    
    %Put into results vector
    plot(Kvec,Counts_norm,'Marker','o','Color',LineColor(ll,:))
 
end

ylabel('P(X=k)');
xlabel('k');
title('Poisson distribution')
legend('\lambda=1','\lambda=4','\lambda=10')
axis([0,20,0,0.4])

% Second method 
% Here is a solution using the PMF: 

% Create a matrix to find the poisson probability of them 
X=0:20;

% The probability of Poisson Distribution for lambda=1
Y1=poisspdf(X,lambda(1));

% The probability of Poisson Distribution for lambda=4
Y2=poisspdf(X,lambda(2));

% The probability of Poisson Distribution for lambda=10
Y3=poisspdf(X,lambda(3));

% Show the result with special color looks like the given link
figure
plot(X,Y1,'Color',[1 0.5 .2],'Marker','o')
hold on
plot(X,Y2,'Color',[0.5 0 .5],'Marker','o')
plot(X,Y3,'Color',[0 0 1],'Marker','o')
hold off
ylabel('P(X=k)');
xlabel('k');
title('Poisson distribution')
legend('\lambda=1','\lambda=4','\lambda=10')
axis([0,20,0,0.4])
