% Luentodemo luennon L6A lopussa (ks. myös luentovideo).
%
% Simulaatio: Kahvikone, jossa hypoteesi H0 on totta / epätotta.
% Miten käy testissä?
% Voit muutella arvoja ja kokeilla mitä tapahtuu.

%#ok<*NOPTS> 

% Hypoteesi
mu0 = 10.0          % Kokeile muuttaa tätä! Entä jos hypoteesi on väärä?

% Merkitsevyystaso
alpha = 0.05        % Kokeile muuttaa tätä!

% Generoiva malli
mu = 10.0           % Kokeile muuttaa tätä!
sigma = 0.5         % Kahvikupin kahvimäärän keskihajonta
n = 30              % Otoksen koko (kuppien lukumäärä)

% Kokeillaan mitä tapahtuu, kun testejä tehdään monta
% (joka kerta n satunnaista kahvikuppia).
%
ntoisto = 1         % Tehtyjen testien määrä. Kokeile muuttaa tätä!
tlist   = [];       % Testisuureen arvot eri kerroilla
hylkylist = logical([]);    % Indikaattorit: Hylättiinkö?
for i=1:ntoisto
    % Uusi data
    x = normrnd(mu, sigma, 1, n);
    
    % Kahvikuppien KESKIARVON keskihajonta.
    % Huom. paljon pienempi kuin kahvikupin keskihajonta.
    std_of_mean = std(x) / sqrt(n);

    % Kupin keskihajontana käytettiin havainnoista laskettua
    % keskihajontaa std(x), joka yleensä poikkeaa sigmasta
    % jonkin verran.  Tarkkaan ottaen pitäisi nyt käyttää
    % t-jakaumaa.

    % Testisuure
    t = (mean(x)-mu0) / std_of_mean;
    tlist = [tlist t];
    
    % P-arvo
    p = 2 * (1-normcdf(abs(t)));
    
    hylky = (p < alpha);
    hylkylist = [hylkylist hylky];
    if hylky
        teksti = 'HYLÄTÄÄN';
    else
        teksti = 'hyväksytään';
    end
    fprintf('Koe %4d: mean=%6.3f, t=%+6.3f, p=%.4f, %s\n', ...
        i, mean(x), t, p, teksti);
end

vasenraja = norminv(alpha/2);
oikearaja = norminv(1 - alpha/2);

clf

xmin = min(-4, min(tlist));
xmax = max(4, max(tlist));
tvals = linspace(xmin, xmax, 1000);
fvals = normpdf(tvals);
plot(tvals, fvals, 'b-', 'linewidth', 2)
hold on
xlabel('t-suureen arvo')
ylabel('t-suureen tiheys jos H0 totta')
for raja=[vasenraja oikearaja]
    plot([raja raja], [0 normpdf(raja)], 'r-', 'linewidth', 2);
end

nhylky = sum(hylkylist);
nhyvak = sum(not(hylkylist));
if length(hylkylist) > 10
    markersize = 5;
else
    markersize = 10;
end
plot(tlist(hylkylist), zeros(1, nhylky), 'ro', 'markerfacecolor', 'r', 'markersize', markersize);
plot(tlist(not(hylkylist)), zeros(1, nhyvak), 'gs', 'markerfacecolor', 'g', 'markersize', markersize);

fprintf('Hylättiin    %3d/%3d kertaa (%6.2f %%)\n', ...
    sum(hylkylist), length(hylkylist), 100*mean(hylkylist));
fprintf('Hyväksyttiin %3d/%3d kertaa (%6.2f %%)\n', ...
    ntoisto-sum(hylkylist), length(hylkylist), 100*mean(1-hylkylist));