% Lecture 6A: Hypothesis testing with coffee machine
%
% Try changing the values and see what happens.
%
% Parameters:
%    nrep   How many tests to run (each with new random cups)
%    n      Number of cups in one test
%    mu0    Mean value that we are testing
%    mu     True mean of the generating distribution
%    sigma  True standard deviation of the generating distribution
%    alpha  Significance level (default = 0.05)
%    use_t  Flag: if true, use the t distribution
%
% If only one test is performed, two plots are made, the top one
% showing the individual cups.
%
function L6Acoffeesimu(nrep, n, mu0, mu, sigma, alpha, use_t)

    if nargin<1, nrep  = 1; end
    if nargin<2, n     = 30; end
    if nargin<3, mu0   = 10; end
    if nargin<4, mu    = 10; end
    if nargin<5, sigma = 0.5; end
    if nargin<6, alpha = 0.05; end
    if nargin<7, use_t  = false; end
    
    tlist   = [];       % Values of test statistic, from each test
    rejectlist = logical([]);    % Indicators (1 = test rejected)
    for i=1:nrep
        % Uusi data
        x = normrnd(mu, sigma, 1, n);
        
        % Observed standard deviation of cups (in this test)
        % using Bessel correction
        sdx = std(x);
    
        % From that, standard deviation of the sample mean.
        std_of_mean = sdx / sqrt(n);
    
        % Compute test statistic
        t = (mean(x)-mu0) / std_of_mean;
        tlist = [tlist t];
        
        % Compute p value
        if use_t
            p = 2 * (1-tcdf(abs(t), n-1));
        else
            p = 2 * (1-normcdf(abs(t)));
        end
        
        reject = (p < alpha);
        rejectlist = [rejectlist reject];
        if reject
            resulttext = 'REJECT';
        else
            resulttext = 'accept';
        end
        fprintf('Test %4d: mean=%6.3f, t=%+6.3f, p=%.4f, %s\n', ...
            i, mean(x), t, p, resulttext);
    end
    
    if use_t
        leftlimit  = tinv(alpha/2,     n-1);
        rightlimit = tinv(1 - alpha/2, n-1);
    else
        leftlimit  = norminv(alpha/2);
        rightlimit = norminv(1 - alpha/2);
    end
    
    clf
    
    if nrep==1
        % Show individual cups in upper panel
        subplot(2,1,1)
        xgrid = linspace(min([mu mu0])-4*sigma, max([mu mu0])+4*sigma, 200);
        fgen  = normpdf(xgrid, mu, sigma);
        fhyp  = normpdf(xgrid, mu0, sigma);
        plot(xgrid, fgen, 'k-', 'linewidth', 4);
        hold on
        plot(xgrid, fhyp, 'b-', 'linewidth', 2);
        plot(x, zeros(size(x)), 'k.', 'markersize', 10);
        plot(mean(x), 0, 'go', 'markerfacecolor', 'g', 'markersize', 10);
        legend('generating', 'hypothesis', 'points', 'sample mean')
        title('Individual cups');

        % Move to lower panel for test statistic
        subplot(2,1,2)
    end

    xmin = min(-4, min(tlist));
    xmax = max(4, max(tlist));
    tvals = linspace(xmin, xmax, 1000);
    if use_t
        fvals = tpdf(tvals, n-1);
    else
        fvals = normpdf(tvals);
    end
    plot(tvals, fvals, 'b-', 'linewidth', 2)
    hold on
    xlabel('value of t statistic')
    ylabel('density, if H0 true')
    for limit=[leftlimit rightlimit]
        if use_t
            plot([limit limit], [0 tpdf(limit, n-1)], 'r-', 'linewidth', 2);
        else
            plot([limit limit], [0 normpdf(limit)], 'r-', 'linewidth', 2);
        end
    end
    title('Test stastistic, distribution vs observed')
    nreject = sum(rejectlist);
    naccept = sum(not(rejectlist));
    if length(rejectlist) > 10
        markersize = 5;
    else
        markersize = 10;
    end
    plot(tlist(not(rejectlist)), zeros(1, naccept), 'gs', 'markerfacecolor', 'g', 'markersize', markersize);
    plot(tlist(rejectlist), zeros(1, nreject), 'ro', 'markerfacecolor', 'r', 'markersize', markersize);
    legend('hypothesis', 'limits');

    fprintf('Rejected    %3d/%3d times (%6.2f %%)\n', ...
        nreject, length(rejectlist), 100*mean(rejectlist));
    fprintf('Accepted    %3d/%3d times (%6.2f %%)\n', ...
        naccept, length(rejectlist), 100*mean(1-rejectlist));

end