
close all
clear all

rng('default');

%
% Q5
%

disp('Q5')

N = 1e6;
X = randn(1,N);

f    = X.^2 + X + 1;
fave = sum(f)/N;

fprintf('average   = %f \n',fave);
fprintf('error     = %f \n',2-fave);

var = 7;
sd  = sqrt(var/N);

fprintf('err bound = %f \n',3*sd);

%
% Q6(b)
%

disp(' '); disp('Q6')

M = 1000;
N = 1000;

fbar = - 2/pi^2;
sig2 = 1/(4*pi^2) + 1/6 - fbar^2;

U = rand(N,1000);
f = U.*cos(pi*U);

%
% Q6(c)
%

Y = sort(sum(f,1)/N);
C = ((1:M) - 0.5)/M;

figure
plot(Y,C, fbar+sqrt(sig2/N)*norminv(C),C)
xlabel('Y'); ylabel('CDF')
legend('Monte Carlo results','Central Limit Theorem',4)

%
% Q6(d)
%

N = 1e6;
U = rand(1,N);
f = U.*cos(pi*U);

S1 = cumsum(f);
S2 = cumsum(f.^2);

NN   = 1e3:1e3:1e6;
mean = S1(NN)./NN;
sd   = sqrt( (S2(NN)./NN - mean.^2)./(NN-1) );

figure
plot(NN,mean, NN,mean-3*sd, NN,mean+3*sd, NN,fbar)
xlabel('N'); ylabel('value')
legend('mean','lower','upper','exact')


