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

%
% Simulation der Gleichverteilung auf dem d-Simplex mittels Verwerfungsmethode
%

%% Zunächst in 3D
d = 3;

% Anzahl und Matrix für Realisierungen
M = 5e3;
X = 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;
   x = rand(1,3);
   if ( min(x>=0) && sum(x) <= 1),
       m = m+1;
       X(m,:) = x;
       TryCount(m) = iter;
       iter = 0;
   end
end
    
figure();
plot3(X(:,1),X(:,2),X(:,3),'.');
xlim([0,1]); ylim([0,1]); zlim([0,1]);
axis equal
title(sprintf('%i Realisierungen der Gleichverteilung auf dem Einheitssimplex in R^3',M))
set(gca,'FontSize',14)
mean(TryCount)


%% Jetzt von 2D bis 6D
M = 1e3;

% Vektor für Anzahl von Versuchen bis zum Erfolg
ds = 2:7;
TryCount = zeros(M, length(ds));

for d = ds,
    % Schleife bis Anzahl Realisierungen erreicht
    m = 0;
    iter = 0;
    while m < M,
       iter = iter + 1;
       x = rand(1,d);
       if ( min(x>=0) && sum(x) <= 1),
           m = m+1;
           TryCount(m,d-1) = iter;
           iter = 0;
       end
    end
end    
figure();
semilogy(ds, mean(TryCount,1),'o-k','LineWidth',2)
xlabel('dimension d')
ylabel('Average Number of Trials')
title('Average Number of Trials until Acceptance for Rejection Sampling for d-Simplex')

% Exakte mittlere Anzahl an Fehlversuchen
hold on;
semilogy(ds, factorial(ds),'+--r','LineWidth',2)
legend('experimental','theoretical')
set(gca,'FontSize',14)