%% UQ-UE3, Aufgabe 1
clear all; clc; close all
rng(17)

%
%% Simulation der Cauchy(0,1)-Verteilung
%

% Vektoren für empirischen Mittelwert und Standardabweichung für M = M_j
MC = zeros(1,11);
STDs = zeros(1,11);

M0 = 100;

% Schleife über Anzahl der Realisierungen
for j = 0:25,
    M = 2^j *M0;
    
    % Simulation mittels Inversionsmethode
    U = rand(M,1);
    X = tan(pi* (U-0.5));
    
    % Monte Carlo Schätzer: empirischer Mittelwert
    MC(1+j) = mean(X);
    
    % empirische Varianz
    STDs(1+j) = std(X);
end

%% Plotten der exakten Dichte und der empirischen relativen Häufigkeiten

% Mittleres 95%-Intervall
int_low = tan(pi * (0.025 - 0.5));
int_upp = tan(pi * (0.975 - 0.5));

% Abschneiden der "Ausreißer"
X_noOutlier = X( find( int_low < X & X < int_upp) );

% absolute Häufigkeiten
nBins = 1000;
[n,x] = hist(X_noOutlier, nBins);

% Plotten der relativen Häufigkeiten gewichtet mit 1/Binbreite 
% (damit Integral/Summe 1)
figure();
bar(x, n/M / (int_upp-int_low) * nBins );

% Plotten der exakten Dichte
hold on;
plot(x, 1/pi * 1./(1+x.^2),'-r','LineWidth',3)
legend('Relative Häufigkeiten', 'Exakte Dichte')
title('Simulation der Cauchy-Verteilung')
set(gca,'FontSize',22)


%% MC-Schätzer für Erwartungswert und Konfidenzintervalle
MC

% Plotten der 95%-Konfidenzintervalle
figure();
for j=0:10,
    M = 2^j * M0;

    semilogx([M,M], [MC(1+j) - 2 * STDs(1+j)/sqrt(M),MC(1+j) + 2 * STDs(1+j)/sqrt(M)], '*-r',...
        M, MC(1+j),'+k','LineWidth',3);
    hold on;
end 
set(gca,'FontSize',22)