function RNG

    clc

    %teil1;

    teil2;

end



function teil1

%% Quasi- und Pseudozufallszahlen

    disp('Mersenne-Twister')

    rand('twister', 5489)

    rand(1,5)

    rand(1,5)

    rand('twister', 5489)

    rand(1,5)



    %Linearer Kongruenzgenerator

    disp('Linearer Kongruenzgenerator')

    rand('seed', 5489)

    rand(1,5)

    rand(1,5)

    rand('seed', 5489)

    rand(1,5)



    %Sobol-Quasi-Zufallszahlen

    disp('Sobol-Quasi-Zufallszahlen')

    p = sobolset(1);

    % nehme die ersten 1000 Werte

    u1 = net(p, 1000);

    u2 = rand(1,1000);

    

    %Plotten der Ergebnisse

    figure('Name','Quasi vs. Pseudo RNs')

    subplot(1,2,1)

    plot(u1,'.')

    title('Sobol-Quasi-Zufallszahlen')

    subplot(1,2,2)

    plot(u2,'.')

    title('Pseudozufallszahlen')



    % Zweidimensionale Sobol-Quasi-Zufallszahlen

    disp('2d-Sobol-Quasi-Zufallszahlen')

    p = sobolset(2);

    v2 = net(p, 1000);

    figure('Name','2d Sobol-Quasi-Zufallszahlen')

    plot(v2,'.')



    

    v2(1,1)

    v2(2,1)



    u2t = v2(:,1);

    diff = u2t - u1;

    disp('  Differenz (erste 10 Werte):')

    diff(1:10)



    v1 = rand(1000,2);

    hold on

    plot(v1,'.r')

    

    

    disp('Weiter mit Enter')

    pause

    

    % Zweidimensionale Normalverteilung

    disp('Zweidimensionale Normalverteilung')

    mu1 = 0;

    mu2 = 0;

    s11 = 10;

    s22 = 10;

    rho = 0.5;

    

    [X,Y] = meshgrid(linspace(-10,10,41));



    function z = f(x1,x2)

        term1  = 1/(2*pi*sqrt(s11*s22*(1-rho^2)));

        term2  = -1/(2*(1-rho^2));



        term3  = (x1-mu1).^2/s11;

        term4  = (x2-mu2).^2/s22;



        term5  = -2*rho*((x1-mu1)*(x2-mu2))/(sqrt(s11)*sqrt(s22));

        z = term1*exp(term2*(term3+term4-term5));

    end

    figure('Name','2d-Normalverteilung')

    surf(X, Y, arrayfun(@f, X, Y))

    title('Zweidimensionale Normalverteilung')

    

    disp('Weiter mit Enter')

    pause



    % Inversmethode

    disp('Inversionsmethode fuer Binomialverteilung')

    figure('Name','Inversionsmethode')

    plot(0:3, binocdf(0:3,3,.5),'O')

    ylim([-0.1 1.1])

    xlim([-0.1 3.1]) 



    p = binocdf(0,3,.5);

    line([0 1],[p p])

    p = binocdf(1,3,.5);

    line([1 2],[p p])

    p = binocdf(2,3,.5);

    line([2 3],[p p])



    %Inverse

    u = 0.65,

    i = binoinv(u, 3, .5);

    line([0 i],[u u],'Color','r')

    line([i i],[u 0],'Color','r','LineStyle',':')



    disp('Weiter mit Enter')

    pause



    %Gewichte einer Poissonverteilung

    disp('Gewichte einer Poissonverteilung')

    bar(0:15, poisspdf(0:15, 3.5), 0.1)

    title('Gewichte der Poissonverteilung')

    xlabel('k')

    ylabel('p_k')

end



function teil2

%% Monte Carlo Methoden

    disp('Monte Carlo Methoden')

    

    

    function y = f(x)

        y = 3*x.^2;

    end

    

    disp('Integration:');

    %Integration mit Simpsonmethode

    quadv(@f,0,1)

    

    function y = mcint(f, a, b, n)

        %Monte Carlo Simulation siehe Script

        ran = a + (b-a).*rand(n,1);

        fran = f(ran);

        y = mean(abs(b-a)*fran);

    end



   % function y = itmci(m)

   %     %Iteration ueber mehrere Simulationen

   %     x = zeros(1,m);

   %     for i=1:m

   %         x(i) = mcint(@f,0,1,1000);

   %     end

   %     y = mean(x);

   % end



    function y =  mcquasiint(f, a,b, n)

        %Monte Carlo mit Quasi-Zufallszahlen

        p = sobolset(1);

        ran = net(p, n) * abs(b-a);

        fran = f(ran);

        y = mean(abs(b-a)*fran);

    end

    

    %runs mit n=100, 1000, 10000

    disp('Aufrufe der Simulationen (mcint vs. mcquasiint)')

    disp('n=100');

    mcint(@f,0,1,100)

    mcquasiint(@f,0,1,100)

    disp('n=1000');

    mcint(@f,0,1,1000)

    mcquasiint(@f,0,1,1000)

    disp('n=10000');

    mcint(@f,0,1,10000)

    mcquasiint(@f,0,1,10000)

    

    disp('mean/std für mcint');

    iterate = zeros(1,100);

    for i=1:100

        iterate(i) = mcint(@f,0,1,10000);

    end

    mean(iterate)

    std(iterate)

    

    disp('mean/std für mcquasiint');

    iterate = zeros(1,100);

    for i=1:100

        iterate(i) = mcquasiint(@f,0,1,10000);

    end

    mean(iterate)

    std(iterate)

    

    

    disp('Weiter mit Enter')

    pause

    

    %Monte Carlo Baender für Min/Max 

    disp('Monte Carlo Baender (Min/Max)')

    m = 100;

    u = 2*rand(m,1);   % Gleichverteilung auf Intervall [0,2]

    g = (1:m)/m;



    stairs(sort(u),g)

    hold on

    title('emp. Verteilungsfunktion mit MC-Baendern')

    ylabel('F')

    

    N = 50;



    U = 2*rand(m, N);

    U = sort(U);

    MCB = [min(U'); max(U')];

    

    stairs(MCB(1,:),g, 'Color','r') 

    stairs(MCB(2,:),g, 'Color','r')

    

    plot(2*g,g,'-.g')    % theor. VF zum Vergleich

    

    % legend('emp VF','MCB','MCB', 'theor VF') 

    hold off

    

 

    

end