%% SS16 UQ UE2, Aufgabe 05

clear all; clc;
rng(17)

%
%% (a) Mittlere Fläche eines zufälliges Dreiecks im Einheitsquadrat
%

% Anzahl an Samples
M = 1e3;    

% Vektor mit Flächen der zufälligen Dreiecke
As = zeros(M,1);

% Schleife zur Simulation und Flächenberechnung der Dreiecke
for m = 1:M,
    % Simulation der zufälligen Eckpunkte
    x1 = rand(1,2); 
    x2 = rand(1,2); 
    x3 = rand(1,2);

    % Flächenberechnung mittels Kreuzprodukt
    As(m) =  0.5 * norm(cross([x2-x1, 0], [x3-x1,0]));
end

% Monte Carlo Schätzer
MC = mean(As);
fprintf('\n MC-Schaezter: %.4f\n', MC);

% geschätzte, asymptotische Konfidenzintervalle
confi_low = MC - norminv(0.975) * std(As)/sqrt(M);
confi_upp = MC + norminv(0.975) * std(As)/sqrt(M);

fprintf('\n asymp. Konfidenzinterval: [%.4f, %.4f] \n', confi_low, confi_upp);

% hierbei:
% 
% std(x)                ...     Standardabweichung des Samples x
% norminv(alpha)     ...     alpha-Quantil der Standardnormalverteilung

% Konfidenzintervalle mittels Chebychev-Ungleichung

cheby_low = MC - 1/sqrt(0.05) * std(As)/sqrt(M);
cheby_upp = MC + 1/sqrt(0.05) * std(As)/sqrt(M);
  
fprintf('\n Tscheby. Konfidenzinterval: [%.4f, %.4f] \n', cheby_low, cheby_upp);

% Konfidenzintervalle mittels Hoeffding-Ungleichung
minA = 0; maxA = 0.5;
hoeff_low = MC - (maxA-minA)/sqrt(M) * sqrt(-log(0.025))/sqrt(2);
hoeff_upp = MC + (maxA-minA)/sqrt(M) * sqrt(-log(0.025))/sqrt(2);
  mean(ok_cheby)
mean(ok_hoeff)

fprintf('\n Hoeff. Konfidenzinterval: [%.4f, %.4f] \n', hoeff_low, hoeff_upp);

% Berry-Esseen-Fehler

BEF = 0.7 * mean(abs(As - mean(As)).^3) / ( std(As)^3 * sqrt(M));

fprintf('\n Berry-Esseen-Fehler: %.4f\n', BEF);


%% Zusatz: Güte der Konfidenzintervalle

ok_asymp = zeros(1,1e4);
ok_cheby = zeros(1,1e4);
ok_hoeff = zeros(1,1e4);
for n = 1:1e4,
    M = 1e3;    
    As = zeros(M,1);
    for m = 1:M,
        % Simulation der zufälligen Eckpunkte
        x1 = rand(1,2); 
        x2 = rand(1,2); 
        x3 = rand(1,2);

        % Flächenberechnung mittels Kreuzprodukt
        As(m) =  0.5 * norm(cross([x2-x1, 0], [x3-x1,0]));
    end

    MC = mean(As);
    
    err = abs(11/144 - MC);
    
    ok_asymp(n) = err <= norminv(0.975) * std(As)/sqrt(M);
    ok_cheby(n) = err <= 1/sqrt(0.05) * std(As)/sqrt(M);
    ok_hoeff(n) = err <= (maxA-minA)/sqrt(M) * sqrt(-log(0.025))/sqrt(2);
end
fprintf('\n Echtes Konfidenzinterval zum asymp. Intervall: %.4f \n', mean(ok_asymp));
fprintf('\n Echtes Konfidenzinterval zum Tscheby.-Intervall: %.4f \n', mean(ok_cheby));
fprintf('\n Echtes Konfidenzinterval zum Hoeffding-Intervall: %.4f \n', mean(ok_hoeff));


%
%% (b) Mittlere Fläche eines zufälliges Dreiecks im Einheitsquadrat,
%     Konfidenzintervalle für verschiedene M
%

% Vektoren für die unteren und oberen Grenzen der Konfidenzintervalle
confi_low = zeros(1,11);
confi_upp = zeros(1,11);

% Äußere Schleife über die Sampleanzahl
M0 = 100;

figure();
for j = 0:10,
    M = 2^j * M0;

    % Vektor mit Flächen der zufälligen Dreiecke
    A = zeros(M,1);

    % Schleife zur Simulation und Flächenberechnung der Dreiecke
    for m = 1:M,
        % Simulation der zufälligen Eckpunkte
        x1 = rand(1,2); 
        x2 = rand(1,2); 
        x3 = rand(1,2);

        % Flächenberechnunorminv(0Size',22).975) * std(As)/sqrt(M)ng mittels Kreuzprodukt
        A(m) =  0.5 * norm(cross([x2-x1, 0], [x3-x1,0]));
    end

    % Monte Carlo Schätzer
    MCS = mean(A);

    % geschätzte, asymptotische Konfidenzintervalle
    confi_low(1+j) = MCS - 2 * std(A)/sqrt(M);
    confi_upp(1+j) = MCS + 2 * std(A)/sqrt(M);
    
    % Zeichnen der Konfidenzintervalle
    plot(j, MCS, 'sg','LineWidth',4); hold on;
    plot([j,j], [confi_low(1+j), confi_upp(1+j)], 's--r','LineWidth',4);

    % Einzeichnen des exakten Wertes
    plot(j, 11/144, 'ko','LineWidth',4);
end
hold off;
xlim([-1,11]); xlabel('j')
legend('MC estimate', 'confi region','truth')
set(gca,'FontSize',22)
%
%% (c) Simulation eines zufälligen Dreiecks auf dem Einheitskreis
%

% Sampleanzahl
M = 1e5;

% zufälliger, gleichverteilter Winkel theta
thetas = 2*pi*rand(M,1);

% Punkte auf dem Kreis
X =  [cos(thetas), sin(thetas)];

% Plotten der Punkte
figure(); subplot(2,2,1);
theta_disk = [0:0.01:2*pi, 0];
plot(cos(theta_disk), sin(theta_disk), '-k', ...
    X(:,1),X(:maxA,2),'.r')
axis equal

% Plotten 3 zufälliger Dreiecke
for i=2:4,
    
    x1 = X(1+3*(i-2),:);
    x2 = X(2+3*(i-2),:);
    x3 = X(3+3*(i-2),:);
    
    subplot(2,2,i);
    plot(cos(theta_disk), sin(theta_disk), '-k', ...
        [x1(1),x2(1),x3(1),x1(1)],[x1(2),x2(2),x3(2),x1(2)],'-r')
    axis equal
end
    
% Plotten der relativen Häufigkeiten und der exakten Verteilungsdichte für 
% x- und y-Koordinate

% absolute Häufigkeiten der x- und y-Koordinate
nBins = 100;
[nx,x] = hist(X(:,1), nBins);
[ny,y] = hist(X(:,2), nBins);

% Zeichnen der relativen Häufigkeiten gewichtet mit 1/Binbreite 
% (damit Integral/Summe 1)
figure();
subplot(1,2,1);
bar(x, nx/M / 2 * nBins ); hold on;

% exakte Dichte
plot(x, 1/pi * 1./sqrt(1-x.^2),'-r','LineWidth',2);
title('Verteilung der x-Koordinate')
xlabel('x')
legend('Relative Häufigkeiten', '1/(\pi * (1-x^2)^{0.5})')

% Analog für y-Koordinate
subplot(1,2,2);
bar(y, ny/M / 2 * nBins ); hold on;
plot(y, 1/pi * 1./sqrt(1-y.^2),'-r','LineWidth',2); % exakte Dichte
title('Verteilung der y-Koordinate')
xlabel('y')
legend('Relative Häufigkeiten', '1/(\pi * (1-y^2)^{0.5})')

% Fläche von M zufaelligen Dreiecken auf dem Einheitskreis
for m = 1:M,
    % Simulation der zufälligen Eckpunkte
    thetas = 2*pi*rand(3,1);
    X =  [cos(thetas), sin(thetas)];

    % Flächenberechnung mittels Kreuzprodukt
    As(m) =  0.5 * norm(cross([X(2,:)-X(1,:), 0], [X(3,:)-X(1,:),0]));
end

% MC-Schaetzung und 95%-Test für Vermutung mu = 3/(2*pi)
MC = mean(As)
minA = 0; 
maxA = 0.75*sqrt(3); % due to Jan's computation
abs(MC - 3/(2*pi)) <= (maxA-minA)/sqrt(M) * sqrt(-log(0.025))/sqrt(2)
