% Sauer Computer Problem 9.1.3 (b)
%
% Gör Monte Carlo typ 1 (uppskatta ett väntevärde). Första delen
% är en illustration av metoden, sedan testar vi den för olika
% antal framslumpade punkter.
%
% Joar Bagge 2016-12-05

%%
% Illustration av Monte Carlo typ 1

clear

P1 = @(x) x.^2 - x + 1/2;
P2 = @(x) -x.^2 + x + 1/2;

figure(1), clf
fplot(P1, [0,1]), hold on
fplot(P2, [0,1])

n = 100; % antal punkter som slumpas fram
X = rand(1,n); % slumpade punkter mellan 0 och 1
% Rita upp punkterna och de längder som vi tar medelvärdet av
for j=1:n
  scatter(X(j), 1/2, 'k')
  plot([X(j), X(j)], [P1(X(j)), P2(X(j))], 'k')
end

% Uppskattningen av väntevärdet är medelvärdet av längderna
A = (1/n) * sum(P2(X)-P1(X))

%%
% Test för flera olika n

ns = 10.^(2:6);
nmax = max(ns);
Xall = rand(1,nmax);
As = [];
for n=ns
  X = Xall(1:n);
  A = (1/n) * sum(P2(X)-P1(X));
  As = [As, A];
end

exact = 1/3; % exakt värde på arean
figure(2)
semilogx(ns, As, '*-') % uppskattningar
hold on
semilogx(ns, 0*ns+exact, 'k--') % exakt svar
title('Integraluppskattning')

figure(3)
loglog(ns, abs(As-exact), '*-') % fel i uppskattningarna
hold on
loglog(ns, 1./sqrt(ns), 'k--') % referenslinje
title('Fel')
% Kör koden flera gånger så blir det fler kurvor i plottarna!
