%% cv06.m 
% A6M33SSL - Statistika a spolehlivost v lekarstvi
% Cviceni 6 - Testovani hypotez (zakladni testy pro jedno rozdeleni +
% parovy pokus)
% 
% 
% 25-10-11, Michal Uricar

clc;
close all; clearvars;

%% Timestamp

fprintf(1,'\n\nStarted on %s\n\n', datestr(now));

%% Chi-squared

fprintf('Kvantily  chi-kvadrat rozdeleni v MATLABu viz prikaz chi2inv.\n');

alpha = 0.05;
confidence = 1-alpha;
qchi2 = zeros(2, 3);   % degrees of freedom in the first row of this matrix
qchi2(1, :) = 18:20;
qchi2(2, :) = chi2inv(confidence, qchi2(1,:));
disp(qchi2);
% or
icdf('chi2', confidence, qchi2(1,:))

%% Chi-squared

for i = 1 : 10
    x = -20:0.1:20;
    y = chi2pdf(x,i);
    z = normpdf(x,0,1);
    figure(1);
    plot(x,y,'-',x,z,'-.')
    legend('Ch-squared distribution', 'normalized normal distribution');
    pause(0.4);
end;

%% Student's t-distribution quantiles

fprintf('Kvantily studentova t-rozdeleni v MATLABu viz prikaz tinv.\n');

alpha = 0.05;
confidence = 1-alpha;
qt = zeros(2, 3);
qt(1, :) = 18:20;
qt(2, :) = tinv(confidence, qt(1,:));
disp(qt);
% or
icdf('t', confidence, qt(1, :))

%%  Comparison of Student's t-distribution and Normalized normal
%%  distribution

for i = 1 : 30
    x = -5:0.1:5;
    y = tpdf(x,i);
    z = normpdf(x,0,1);
    figure(1);
    plot(x,y,'-',x,z,'-.')
    legend('Student''s t-distribution', 'normalized normal distribution');
    title(sprintf('t(%d)', i));
    pause(0.2);
end;

%% Poisson's distribution

fprintf('Kvantily poissonova rozdeleni v MATLABu viz prikaz poissinv.\n');

alpha = 0.05;
confidence = 1-alpha;
qpoiss = zeros(2, 3);
qpoiss(1, :) = 18:20;
qpoiss(2, :) = poissinv(confidence, qpoiss(1,:));
disp(qpoiss);
% or
icdf('poiss', confidence, qpoiss(1, :))

%% Poisson's distribution

for i = 1 : 20
    x = 0:0.001:50;
    y = poisspdf(x,i);
    z = normpdf(x,i,sqrt(i));
    figure(1);
    plot(x,y,'-',x,z,'-.');
    ylim([0, 0.2]);
    legend('Poisson''s distribution', 'normalized normal distribution');
    title(sprintf('lambda = %d', i));
    pause(0.5);
end;

%% Cviceni 1 - Test parametru p rozdeleni Bi(1, p)

alpha = 0.05;
n = 2000;
b = 835;
x = b/n;

T = (x - 0.5)/(sqrt((x*(1-x))/2000));
if (abs(T) > norminv(1-alpha/2, 0, 1))
    fprintf('Zamitame.\n');
else 
    fprintf('Nezamitame.\n');
end;
% dosazena vyznamnost
vyznamnost = 2 * min(normcdf(T, 0, 1), 1-normcdf(T, 0, 1))

%% Cviceni 2 - 12.2.4

n = 100000;
v = 16;
q = v/n;

alpha = 0.05; 

% riziko se nezvysuje => Binomicke rozdeleni Bi(n, q), CLV -> ~ N(nq,nq(1-1))
kriticka_hodnota = norminv(1-alpha, 0, 1);
vzorek = kriticka_hodnota^2*(1-q)/q;

2*(vzorek/n)*v

% overeni
vzorek = round(vzorek);
flag = true;
krit = inf
sum = 0;
k = 0;
while (krit > alpha)
    flag = false;
    sum = sum + nchoosek(vzorek, k)*q^k*(1-q)^(vzorek-k);
    krit = 1 - sum;
    k = k + 1;
end;
k = k-1

%% Cviceni 3 - 12.2.3

x = 5.3;
u = 5;
sx = 0.3;
n = 10;
alpha = 0.01;
T = sqrt(n)*(x-u)/(sx)
if (abs(T) > tinv(1-alpha/2, n-1))
    fprintf('Zamitame.\n');
else
    fprintf('Nezamitame.\n');
end;
% dosazena vyznamnost
vyznamnost = 2 * min(tcdf(T, n-1), 1-tcdf(T, n-1))

%% Cviceni 4 - 12.2.30

teploty = [37.5, 38.2, 38.3, 38.6, 38.2, 37.6, 37.2, 36.9, 36.6;
           37.7, 38.3, 38.5, 38.6, 38.3, 37.7, 37.4, 36.7, 36.5 ];
alpha = 0.05;

delta = teploty(1,:)-teploty(2,:);
n = size(teploty, 2);
dmean = mean(delta);
dvar = var(delta);
dstd = std(delta);  % nebo dstd = sqrt(dvar);
t = sqrt(n)*dmean/dstd

if (abs(t) > tinv(1-alpha/2, n-1))
    fprintf('Zamitame.\n');
else
    fprintf('Nezamitame.\n');
end;

%% Timestamp

fprintf(1,'\n\nFinished on %s\n\n', datestr(now));