% Poni¿sza tabela prezentuje wyniki 9 pacjentów wykonuj±cych pewien test diagnostyczny
%  przed podaniem leku i po podaniu leku. 

%przed	po
A=[1.83	0.878
0.50	0.647
1.62	0.598
2.48	2.05
1.68	1.06
1.88	1.29
1.55	1.06
3.06	3.14
1.3	1.29];
%Skonstruowaæ test, który pozwoli stwierdziæ czy lek jest skuteczny. 
%Porównaæ dwie wersje testu bootstrapow± (losowanie z powtórzeniami) i permutacyjn±.

% Jako miare tego czy lek jest skuteczny przyjmiemy ró¿nicê (po-przed)
% Bêdziemy wierzyli ¿e lek dzia³a jeœli ta ró¿nica jest istotnie ró¿na od
% zera

r=A(:,2)-A(:,1);
disp('srednia roznica:')
mr=mean(r);
disp(mr)
% H0: wszystko jedno ktory pomiar jest przed a ktory po (r_srednie równe
% zero)
% H1: r_srednie ró¿ne zero


% Musimy wytwarzac symulowane r zgodnie z hipoteza zerowa, ze nie ma
% znaczenia która wartosc jest przed, a która po. trzeba tylko uwarzac ¿eby
% nie pomieszac pacjentow

% wersja bootstrpowa: bierzemy pacjentow z powtorzeniami - zak³adamy
% reprezentatywna grupa bardzo duzej populacji

N_rep=1e4;
r_boot=zeros(1,N_rep);
N=length(r);
for i=1:N_rep
    ix=randsample(N,N,true); % wybieramy pacjentow z powtorzeniami
    B=A(ix,:);
    for j=1:N   % mieszamy losowo przypisujemy przypadki do grupy przed i po
        s=rand(1,1);
        if s>0.5
            przed(j)=B(j,1);
            po(j)=B(j,2);
        else
            przed(j)=B(j,2);
            po(j)=B(j,1);
        end
    end
    rr=po-przed;     
    r_boot(i)=mean(rr);
end
subplot(221)
hist(r_boot,30)
disp('dla repróbkowanego testu dwustronnego: ')
p_h0=length(find(abs(r_boot)>abs(mr)))/N_rep;
disp(sprintf('p_H0: %.3f',p_h0))
title(sprintf('test repróbkowany p_{H0}: %.3f',p_h0))


% wersja permutacyjna: korzystamy za ka¿dym razem ze wszystkich dostepnych
% pacjentow

N_rep=1e4;
r_boot=zeros(1,N_rep);
N=length(r);
for i=1:N_rep
    for j=1:N   % mieszamy losowo przypisujemy przypadki do grupy przed i po
        s=rand(1,1);
        if s>0.5
            przed(j)=A(j,1);
            po(j)=A(j,2);
        else
            przed(j)=A(j,2);
            po(j)=A(j,1);
        end
    end
    rr=po-przed;     
    r_boot(i)=mean(rr);
end
subplot(222)
hist(r_boot,30)
disp('dla permutacyjnego testu dwustronnego: ')
p_h0=length(find(abs(r_boot)>abs(mr)))/N_rep;
disp(sprintf('p_H0: %.3f',p_h0))
title(sprintf('test permutacyjny p_{H0}: %.3f',p_h0))

% jesli mielibysmy jakies przeslanki, zeby zalozyc normalnosc naszych
% danych wejsciowych to mozna zastosowac test t 
disp('dla parowanego testu t dwustronnego: ')
subplot(223)
normplot(r)
[h, p]=ttest(A(:,2),A(:,1))
title(sprintf('parowany test t p_{H0}: %.3f',p))

subplot(224)
axis
axis off
disp('Testy nieparametryczne testuja czy mediana roznicy jest 0')
p_st = signtest(A(:,2),A(:,1)); % zal: dane pochodza z dowolnej ciaglej dystrybucji
s1=sprintf(' test znkow p_{H0}: %.3f',p_st)
p_sr = signrank(A(:,2),A(:,1)); % zal: dane pochodza z dowolnej ciaglej symetrycznej dystrybucji 
s2=sprintf(' test rang p_{H0}: %.3f',p_sr)
str(1)={'Testy nieparametryczne testuja'};
str(2)={'czy mediana roznicy jest 0'};
str(3)={' '};
str(4)={s1};
str(5)={s2};
text(0.1, 0.7,str)


