%% UQ-UE3, Aufgabe 3
clear all; clc;

%
% Simulation der Gleichverteilung auf der Kugeloberfläche
%


%
%% (a) Simulation mittels Verwerfungsmethode
%

% Anzahl an Samples
M = 5e3;    

% Vektor mit Realisierungen
S = zeros(M,3);

% Vektor für Anzahl von Versuchen bis zum Erfolg
TryCount = zeros(M,1);

% Schleife bis Anzahl Realisierungen erreicht
m = 0;
iter = 0;
while m < M,
    iter = iter + 1;
 
    % Gleichverteilung auf [-1,1]^3
    x = 2*rand(1,3) - 1; 
    
    % Testen, ob Punkt in Einheitskugel
    if (norm(x) <= 1 )
        m = m+1;
        S(m,:) = x/norm(x);
        TryCount(m) = iter;
        iter = 0;
    end
end

AverageTrialNumber = mean(TryCount)

% Theoretische Akzeptanzwahrscheinlichkeit
p = 4/3 * pi /(2^3);
ExactTrialNumber = 1/p

% Plotten der Realisierungen
sampleplot = figure();
subplot(1,3,1);
plot3(S(:,1),S(:,2),S(:,3),'.k')
title('Mit Verwerfungsmethode')
axis equal

% Plotten der Histogramme für x-,y- und z-Koordinaten
histplot = figure();

subplot(3,3,1); hist(S(:,1),25)
xlabel('x_1')
title('Mit Verwerfungsmethode')

subplot(3,3,2); hist(S(:,2),25)
xlabel('x_2')
title('Mit Verwerfungsmethode')

subplot(3,3,3); hist(S(:,3),25)
xlabel('x_3')
title('Mit Verwerfungsmethode')


%
%% (b) Direkte Simulation mittels multivariater Standardnormalverteilung
%

% Anzahl an Samples
M = 5e3;    

% Vektor mit Realisierungen
S = randn(M,3);
S = S ./ repmat( sqrt(sum(S.^2,2)),1,3); 

% Plotten der Realisierungen
figure(sampleplot);
subplot(1,3,2);
plot3(S(:,1),S(:,2),S(:,3),'.k')
title('Mit direkter Methode')
axis equal

% Plotten der Histogramme für x-,y- und z-Koordinaten
figure(histplot);
subplot(3,3,4); hist(S(:,1),25)
xlabel('x_1')
title('Mit direkter Methode')

subplot(3,3,5); hist(S(:,2),25)
xlabel('x_2')
title('Mit direkter Methode')

subplot(3,3,6); hist(S(:,3),25)
xlabel('x_3')
title('Mit direkter Methode')


%
%% (c) Untersuchung alternativer, direkter Simulation
%

S2 = zeros(M,3);

for m = 1:M,
    % Gleichverteilung auf [-1,1]^3
    u1 = 2 * rand(1,1) -1;
    u2 = 2 * rand(1,1) -1;
    u3 = 2 * rand(1,1) -1;

    x = [u1, u2, u3];

    S2(m,:) = x/norm(x);
end

figure(sampleplot);
subplot(1,3,3);
plot3(S2(:,1),S2(:,2),S2(:,3),'.k')
title('Alternative Methode')

axis equal

figure(histplot);
subplot(3,3,7); hist(S2(:,1),25)
xlabel('x_1')
title('Alternative Methode')

subplot(3,3,8); hist(S2(:,2),25)
xlabel('x_2')
title('Alternative Methode')

subplot(3,3,9); hist(S2(:,3),25)
xlabel('x_3')
title('Alternative Methode')

%
%% Zusatz: Verteilung der Koordinaten für Gleichverteilung auf Einheitssphäre im R^d
% (dauert einige Minuten...)

% Anzahl an Samples
M = 5e4;    

% Anzahl an Bins für Plot von relativen Häufigkeiten
nBins = 100;

% Normierungskonstanten für vermutete exakte Dichten:
% f_d(x) ~ (1 - x^2) ^( (d-3)/2 ) 
%
% (provided by WolframAlpha)

cs(2) = pi;
cs(3) = 2;
cs(4) = pi/2;
cs(5) = 4/3;
cs(6) = 3/8 * pi;
cs(7) = 16/15;
cs(8) = 5/16 * pi;
cs(9) = 32/35;
cs(10) = 35/128 * pi;

% Schleife über Dimension
figure();
for d = 2:10,

    % Vektor mit Realisierungen
    S = zeros(M,d);

    % Schleife bis Anzahl Realisierungen erreicht
    m = 0;
    while m < M,
        % Gleichverteilung auf [-1,1]^d
        x = 2*rand(1,d) - 1; 

        % Testen, ob Punkt in Einheitskugel
        if (norm(x) <= 1 )
            m = m+1;
            S(m,:) = x/norm(x);
        end
    end
    % Plotten der relativen Häufigkeiten
    subplot(3,3,d-1);
    [n,x] = hist(S(:,1), nBins);
    bar(x, n/M / 2 * nBins );
    % Plotten der vermuteten dichte
    hold on;
    plot(x, 1/cs(d) * sqrt((1-x.^2)).^(d-3),'-r','LineWidth',2)
    title(sprintf('Verteilung der x_1-Koordinate fuer Uni(S^{%i-1})',d))
end
