Statystyka w R


Idź do treści

Testy parametryczne


Kolejnym (obok teorii estymacji) działem wnioskowania statystycznego jest weryfikacja hipotez statystycznych. Testy statyczne dzielą się na parametryczne (dotyczące nieznanych parametrów populacji) i nieparametryczne. Weryfikacja hipotez w wielu wypadkach wymaga spełniania założenia normalności rozkładu.


####################################################################
############ Testy parametryczne - wybrane zagadnienia #############
####################################################################

### Testy dla jednej próby ###

### Test dla średniej - znane sigma
test_srednia_sig=function(x,m0,sigma,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić, czy spełnione są przyjęte założenia")
n=length(x)
stat_test=(mean(x)-m0)*sqrt(n)/sigma # wartość statystyki testowej
FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}
# Przykład:
dane=rnorm(10,100,15)
test_srednia_sig(dane,120,15,0.05,H1="W")

### Test dla średniej - nieznane sigma, duża próba
test_srednia_duza=function(x,m0,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić liczebności prób i czy spełnione są przyjęte założenia")
n=length(x)
S2=var(x)*(n-1)/n
S=sqrt(S2)
stat_test=(mean(x)-m0)*sqrt(n)/S
FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}
# Przykład:
dane=rnorm(50,120,15)
test_srednia_duza(dane,120,0.05,H1="W")

### Test dla średniej - nieznane sigma, mała próba
test_srednia_mala=function(x,m0,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić, czy spełnione są przyjęte założenia")
n=length(x)
S2=var(x)*(n-1)/n
S=sqrt(S2)
stat_test=(mean(x)-m0)*sqrt(n-1)/S
tW_alfa=qt(alfa, (n-1), lower.tail = FALSE)
tM_alfa=qt(alfa, (n-1), lower.tail = TRUE)
tR_alfa=qt((alfa/2), (n-1), lower.tail = FALSE)

if (H1=="W") t_alfa=tW_alfa
if (H1=="M") t_alfa=tM_alfa
if (H1=="R") t_alfa=tR_alfa

if (abs(stat_test)>=abs(t_alfa))wynik="H1"
if (abs(stat_test)<abs(t_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,t_alfa=t_alfa,wynik=wynik)
}
#Przykład:
dane=rnorm(10,120,15)
test_srednia_mala(dane,130,0.05,H1="M")
#powyższy test można również przeprowadzić następująco:
t.test(dane-130,conf.level = 0.95,alternative="less")

### Test dla frakcji
test_frakcja=function(m,n,p0,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
if (n<=100) stop("liczebność próby musi większa od 100")
est=m/n
stat_test=(est-p0)/(sqrt(p0*(1-p0)/n))
FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}

# Przykład 1
palacy=50 #liczba palących w próbie
wszyscy=200 #liczebność próby
test_frakcja(palacy,wszyscy,0.5,0.01,H1="M")
# Przykład 2
internet=45 #liczba posiadających dostęp do internetu w próbie
wszyscy=90 #liczebność próby
test_frakcja(internet,wszyscy,0.6,0.01,H1="R")


#####################################################

### Testy dla dwóch prób ###

### test równości średnich - znane sigmy
test_2srednie_sig=function(x,y,sigma_x,sigma_y,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić, czy spełnione są przyjęte założenia")
n1=length(x)
n2=length(y)
sigma2_1=sigma_x^2
sigma2_2=sigma_y^2
stat_test=(mean(x)-mean(y))/sqrt((sigma2_1/n1)+(sigma2_2/n2))
FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}
# Przykład:
dane_x=rnorm(10,100,15)
dane_y=rnorm(20,120,20)
test_2srednie_sig(dane_x,dane_y,15,20,0.05,H1="R")

### test równości średnich - nieznane sigmy, duże próby
test_2srednie_duza=function(x,y,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić liczebności prób i czy spełnione są przyjęte założenia")
n1=length(x)
n2=length(y)
S2_1=var(x)*(n1-1)/n1
S2_2=var(y)*(n2-1)/n2
stat_test=(mean(x)-mean(y))/sqrt((S2_1/n1)+(S2_2/n2))
FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}

dane_x=rnorm(40,100,15)
dane_y=rnorm(60,120,20)
test_2srednie_duza(dane_x,dane_y,0.05,H1="R")


### test równości średnich - nieznane sigmy, małe próby
test_2srednie_mala=function(x,y,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
zal=c("Sprawdzić liczebności prób i czy spełnione są przyjęte założenia")
n1=length(x)
n2=length(y)
S2_1=var(x)*(n1-1)/n1
S2_2=var(y)*(n2-1)/n2
stat_test_a=mean(x)-mean(y)
stat_test_b=(n1*S2_1+n2*S2_2)/(n1+n2-2)
stat_test_c=stat_test_b*((1/n1)+(1/n2))
stat_test=stat_test_a/sqrt(stat_test_c)

tW_alfa=qt(alfa, (n1+n2-2), lower.tail = FALSE)
tM_alfa=qt(alfa, (n1+n2-2), lower.tail = TRUE)
tR_alfa=qt((alfa/2), (n1+n2-2), lower.tail = FALSE)

if (H1=="W") t_alfa=tW_alfa
if (H1=="M") t_alfa=tM_alfa
if (H1=="R") t_alfa=tR_alfa

if (abs(stat_test)>=abs(t_alfa))wynik="H1"
if (abs(stat_test)<abs(t_alfa))wynik="H0"
list(zal=zal,stat_test=stat_test,t_alfa=t_alfa,wynik=wynik)
}
#Przykład:
dane_x=rnorm(10,100,15)
dane_y=rnorm(20,120,20)
test_2srednie_mala(dane_x,dane_y,0.05,H1="R")
# ten sen problem, ale innym testem:
t.test(dane_x,dane_y,alternative = "two.sided", conf.level = 0.95)


### test równości dwóch frakcji
test_2frakcje=function(m1,n1,m2,n2,alfa,H1) {
#H1 - W, M, R - większe, mniejsze, różne
if ((n1<=100)| (n2<=100)) stop("liczebności prób muszą być większe od 100")
n_=n1*n2/(n1+n2)
p_=(m1+m2)/(n1+n2)
q_=1-p_

stat_test_a=(m1/n1)-(m2/n2)
stat_test_b=p_*q_/n_
stat_test=stat_test_a/sqrt(stat_test_b)

FuR=1-(alfa/2)
uR_alfa=qnorm(FuR)
FuW=1-alfa
uW_alfa=qnorm(FuW)
FuM=alfa
uM_alfa=qnorm(FuM)
if (H1=="W") u_alfa=uW_alfa
if (H1=="M") u_alfa=uM_alfa
if (H1=="R") u_alfa=uR_alfa

if (abs(stat_test)>=abs(u_alfa))wynik="H1"
if (abs(stat_test)<abs(u_alfa))wynik="H0"
list(stat_test=stat_test,u_alfa=u_alfa,wynik=wynik)
}

#Przykład1
kobiety_palace=30 # liczba kobiet palących w próbie
kobiety=120 # liczebność próby wylosowanej z populacji kobiet
mezczyzni_palacy=60 # liczba mężczyzn palących w próbie
mezczyzni=200 # liczebność próby wylosowanej z populacji mężczyzn
test_2frakcje(kobiety_palace,kobiety,mezczyzni_palacy,mezczyzni,0.05,H1="R")

#Przykład2
kobiety_palace=20 # liczba kobiet palących w próbie
kobiety=90 # liczebność próby wylosowanej z populacji kobiet
mezczyzni_palacy=60 # liczba mężczyzn palących w próbie
mezczyzni=200 # liczebność próby wylosowanej z populacji mężczyzn
test_2frakcje(kobiety_palace,kobiety,mezczyzni_palacy,mezczyzni,0.05,H1="R")






Powrót do treści | Wróć do menu głównego