1.Utworzyć na trzy sposoby: (c, :, seq) wektor składający się z elementów 4,3,2,1,0.

w=c(4,3,2,1,0)

v=c(4:0)

u=seq(4:0)

2.Wypisać (możliwie jak najkrótszą instrukcją) trzecie potęgi liczb od 1 do 20.

(1:20)^3

3.Wypisać (możliwie jak najkrótszą instrukcją) kwadraty liczb parzystych od 10 do 30.

seq(10,30,2)^2=09

4.Podać (możliwie jak najkrótszą instrukcją) sumę sześcianów liczb od 1 do 98.

5.(*) Obliczyć średni zarobek dla danych z pliku Wyk_Zarobki.xls

dane=read.csv2("dane2.csv")

w=c(0,seq(500,10500,1000))

u=dane[,4]

S=0

for(n in 1:305){

if(u[n]>1) S=S+(u[n]-2)*1000+500}

cat(S/305)

Dla danych zawartych w pliku wzrost.xls (próbka statystyczna opisująca wzrost dorosłego mężczyzny w Polsce) wyznaczyć podstawowe statystyki opisowe:

* mediana, kwartyle, rozstęp;

* średnia, wariancja - wraz z przedziałem ufności na poziomie 0,9;

* skośność i kurtoza

wzrost=read.csv2("wzrost.csv")

summary(wzrost)

quantile(wzrost[,1], 0.25)

quantile(wzrost[,1], 0.75)

median(wzrost[,1])

sd(wzrost[,1])

(r=max(wzrost)-min(wzrost))

kurtosis(wzrost[,1],type=2) #nieobiciazony bo type2 i e1071

#Skosnosc

n=length(wzrost[,1])

s=n/(n-1)/(n-2)*sum(wzrost[,1]-mean(wzrost[,1])^3)/(sd(wzrost[,1))^3)

s=e1071::skewness(wzrost[,1],type=2) #trzeba przez e1071 obliczyc

#przedział ufności dla wartosci oczekiwanej

a=0,1 #błąd zmienic alfa

qt(a,n-1,lower.tail=T)

(m=mean(wzrost[,1])-qt(1-a/2,n-1,lower.tail=T)*sd(wzrost[,1])/(n)^(1/2))

n=mean(wzrost[,1])+qt(1-a/2,n-1,lower.tail=T)*sd(wzrost[,1])/(n)^(1/2)

cat("przedział ufnośći dla wartosci oczekiwanej wynosi:","[",m,",",n, "]")

#Przedział ufnosci dla wariancji

b=n*(sd(wzrost[,1])^2)/(qchisq(1-a/2,n-1,lower.tail=T))

c=n*(sd(wzrost[,1])^2)/(qchisq(a/2,n-1,lower.tail=T))

cat("przedział ufnośći dla wariancji wynosi:", "[",b,",",c, "]")

skewness(wzrost)

Dla danych zawartych w pliku zakupy.xls (próbka opisująca kwotę zakupów w sklepie wielkopowierzchniowym w zależności od pewnych cech klienta) obliczyć wybrane statystyki z podziałem na kobiety i mężczyzn:

* średnia, wariancja;

* mediana, kwartyle, rozstęp;

* skośność i kurtoza

zakupy=read.csv2("zakupy.csv")

N=length(zakupy[,3])

zakupy[,2]

w=(zakupy[,2]=="K")

zakupy$WYDATEK[w]

kobiety=zakupy$WYDATEK[zakupy[,2]=="K"]

man=zakupy$WYDATEK[zakupy[,2]=="M"]

mean(kobiety)

mean(man)

var(kobiety)

var(man)

median(kobiety)

median(man)

quantile(kobiety,0.25)

quantile(kobiety,0.5)

quantile(kobiety,0.75)

quantile(man,0.25)

quantile(man,0.5)

quantile(man,0.75)

skewness(kobiety,type=2)

skewness(man,type=2)

kurtosis(kobiety,type=2)

kurtosis(man,type=2)

Wyznaczyć wartość oczekiwaną i wariancję rozkładu podanego w tabeli:

wartość

1 2 3 4 5

p-stwo

7/30 1/6 4/15 2/15 1/5

X=c(1,2,3,4,5)

P=c(7/30, 1/6, 4/15, 2/15, 1/5)

wartOczekiwana=function(P,X){

EX=0

for (i in 1:5) {

EX=EX+X[i]*P[i]

}

return(EX)}

wartOczekiwana(P,X)

EXX=function(P,X) {

EXX=0

for(i in 1:5){

EXX=EXX+X[i]^2*P[i]}

return(EXX)}

EXX(P,X)

Wariancja=function(P,X) {

VarX=EXX(P,X)-wartOczekiwana(P,X)^2

return(VarX) }

cat("Wariancje jest równa ",Wariancja(P,X))

Napisać funkcję pu_test(n, mi, sigma, poziom), która:

* losuje n liczb z rozkładu normalnego N(mi, sigma)

* w oparciu o uzyskane dane wyznacza przedział ufności [x, y] dla średniej na poziomie poziom (przyjmując nieznaną wariancję)

* zwraca wartość TRUE, jeśli mi wpada do wyznaczonego przedziału i FALSE w przeciwnym razie.

Utworzoną funkcję wykorzystać do następującego eksperymentu: policzyć ile razy wypadnie fałsz przy 1000 wykonaniach dla n=100 i poziom = 0.9.

pu_test=function(n,u,q,a){

w=c(rnorm(n,mean=u,sd=q))

m=mean(w)-qt(1-(1-a)/2,n-1,lower.tail=T)*((q/n)^(1/2))

n=mean(w)+qt(1-(1-a)/2,n-1,lower.tail=T)*((q/n)^(1/2))

if(u>=m || u<=n)

y=TRUE

else

y=FALSE

y}

pu_test(100,5,34,0.9)

ileliczb=function(n,u,q,a) {

S=0; F=0

for(i in 1:1000) {

if(pu_test(n,u,q,a)==TRUE)

S=S+1

else

F=F+1}

cat(" TRUE otrzymano:",S,"razy\n","FALSE otrzymano:",F,"razy\n")}

ileliczb(100,0,1,0.9)

Napisać funkcję statystyki(x), która dla zadanego wektora x będzie wyznaczać podstawowe statystyki opisowe z jego współrzędnych (przynajmniej średnią, odch. standardowe, skośność, kurtozę, mininum, maksimum). Następnie wylosować po 100 liczb z rozkładu N(0,1) zapisując wyniki w zmiennych x i y oraz wykonać tę funkcję dla x, y, x+y, x-y, p(x)+p(y), p(x)-p(y), gdzie p(w) oznacza wektor zawierający uporządkowane współrzędne wektora w.

statystyki=function(x) {

n=length(x)

średnia=mean(x)

odchylenieStand=sd(x)

min=min(x)

max=max(x)

kurtoza=e1071::kurtosis(x)

skośność=e1071::skewness(x)

cat("Średnia:",średnia,"\n")

cat("Odchylenie standardowe:",odchylenieStand,"\n")

cat("Minimum:",min,"\n")

cat("Maksimum:",max,"\n")

cat("Kurtoza:",kurtoza,"\n")

cat("Skośność:",skośność,"\n")}

x=c(2,3,4,5,1,2,2,2,2)

statystyki(x) #czym sie roznia funkcje return od cat

x=c(rnorm(100,mean=0,sd=1))

y=c(rnorm(100,mean=0,sd=1))

statystyki(x)

statystyki(y)

statystyki(x+y)

#nierosnaca dodajemy za rpzecinkiem slowo TRUE

statystyki(sort(x)+sort(y))

statystyki(sort(x)-sort(y))

Napisać funkcję srednia_p(x, p) obliczającą średnią potęgową rzędu p elementów wektora x. Następnie dla 1500 liczb pseudolosowych, wylosowanych z rozkładu jednostajnego na przedziale [1,10], obliczyć średnią arytmetyczną, średnią harmoniczną, średnią geometryczną oraz średnią kwadratową (potęgową rzędu 2).

średnia_p=function(x,p)

{sum=0

n=length(x)

for(i in 1:n) {

sum=sum+((x[i])^p)}

cat("średnia rzędu",p,"wynosi:", (sum/n)^(1/p),"\n")}

x=c(2,3,4,5,1,2,2,2,2)

średnia_p(x,3)

średnia_p(runif(1500,1,10),1)

średnia_p(runif(1500,1,10),-1)

średnia_p(runif(1500,1,10),2)

y=c(runif(1500,1,10))

M1=prod(y[1:300])^(1/length(y))

M2=prod(y[301:600])^(1/length(y))

M3=prod(y[601:900])^(1/length(y))

M4=prod(y[901:1200])^(1/length(y))

M5=prod(y[1201:1500])^(1/length(y))

cat("średnia geometryczna wynosi:",prod(M1,M2,M3,M4,M5))

#prościej:

cat("średnia geometryczna wynosi:",prod(y^(1/length(y))))

#jeszcze inny sposob- tak liczy R:

exp(mean(log(y)))

Kierowca Formuły 1 podczas wyścigu o Grand Prix Monako (na torze Monte Carlo, a zresztą, czy ma to jakieś znaczenie?) przejeżdżał kolejne okrążenia toru ze średnimi prędkościami podanymi w pliku Formula1.csv. Z jaką średnią prędkością przejechał cały wyścig?

formula=read.csv2("Formula1.csv",header=F)

statystyki(formula)

z=c(formula)

średnia_p(formula,-1)

$moja formula poprawiona

średnia_p=function(x,p){

sum=0

n=length(x)

for(i in 1:n) {

sum=sum+((x[i,])^p)}

cat("średnia rzędu",p,"wynosi:", (sum/n)^(1/p),"\n")}

średnia_p(c(formula),-1)

srednia_p=function(x, p){

a=(sum(x^p)/length(x))^(1/p)

cat(a)}

Sporządzić wykres ramkowy (pudełkowy) dla wieku nowożeńców (dane Sluby.csv).

(*) Czy liczba 25 wpada do przedziału ufności (na poziomie 95%) dla mediany wieku kobiety? A mężczyzny?

sluby=read.csv2("sluby.csv") #(a)

boxplot(sluby); #(b)

i=25; #badana liczba

alfa=0.05; #poziom istotnosci

N=length(sluby$ZONA); #liczebnosc probki

q=0.5; #mediana

a=ceiling(N*q-qnorm(1-alfa/2)*sqrt(N*q*(1-q)));

b=ceiling(N*q+qnorm(1-alfa/2)*sqrt(N*q*(1-q))); #przedzial ufnosci dla mediany [a,b]

test=function(){

if (sort(sluby$ZONA)[a]<=i && i<=sort(sluby$ZONA)[b])

cat("Liczba",i,"wpada do przedzia=B3u ufnosci dla zon\n");

if (sort(sluby$MAZ)[a]<=i && i<=sort(sluby$MAZ)[b])

cat("Liczba",i,"wpada do przedzialu ufnosci dla mezow\n");}

test();

#liczba 25 z prawdopodobienstwem okolo 0,95 stoi na pozycji miedzy =a i b,=20

gdzie mediana dla zon jest 23, dla mezow 24

Dla danych z pliku wzrost.csv sporządzić wykres zawierający:

* histogram (około 15 klas)

* jądrowy estymator gęstości

* gęstość rozkładu normalnego o parametrach estymowanych na podstawie danych

wzrost=read.csv("wzrost.csv",header=3DF)

hist(wzrost[,1],15,prob=T,col="blue",main="Histogram wzrostu") =

#prob=T zmieniam jednostke...

lines(density(wzrost[,1]),col="red",lwd=3)

#bez prob=t nie narysuje dwoch funkcji....

plot(density(wzrost[,1]),main="Jadrowy estymator gestosci")

Dane z pliku wzrost.csv rozłożyć w szereg rozdzielczy:

* o 15 klasach taki, że minimum jest środkiem pierwszego, a maksimum środkiem ostatniego przedziału

* (*) o 6 klasach, różnej długości, za to, w przybliżeniu, jednakowej liczności. W oparciu o ten szereg narysować histogram.

k=15;

min(wzrost[,1])

dl=3D(max(wzrost[,1])-min(wzrost[,1]))/(k-1)

podzial=dl*0:k-dl/2+min(wzrost[,1])

table(cut(wzrost[,1],podzial))

hist(wzrost[,1],podzial,col="blue",main="Histogram wzrostu") #(b)

k=6;

d=signif(length(wzrost[,1])/(k),6)

sort(wzrost[,1])[1+d*0:k] #poczatki/krance odcinkow podzialu

podzial2=sort(wzrost[,1])[1+d*0:k];

hist(wzrost[,1],podzial2,col="yellow",main="Histogram wzrostu")

przyklad do wykładu 5

wzrost = read.csv2("Wzrost.csv", header = F)[,1]

?t.test

?chisq.test

apropos("test")

#testuje, czy przecietny wzrost jest równy 170, przeciwko zaprzeczeniu

t.test(wzrost, mu = 170)

#No dobrze, ale czy wzrost ma rozkład normalny?

#test chi^2 - zgodności

mi = mean(wzrost)

sigma = sd(wzrost)

ile = length(wzrost)

#15 klas, potem 20

k = 15

w = 0:k

#1 podejscie - jednakowe p-stwa

podzial = qnorm(w/k, mi, sigma)

podzial

podzial[1] = min(wzrost)-5

podzial[length(podzial)] = max(wzrost)+5

szereg = table(cut(wzrost, podzial))

szereg

(t=chisq.test(szereg))

#jak sie dobrać do wyników?

names(t)

#ale estymowałem parametry (dwa) rozkładu normalnego

#poprawiona p-wartość

1-pchisq(t$statistic, k-3)

print(hist(wzrost, podzial, col = "yellow"))

szereg

#oczywiscie moge sam policzyć:

teor=ile/k

(chi2=sum((szereg-teor)^2/teor))

#jak wiec sobie poradzić?

#może jednak podzial od strony danych, 20 klas

k = 20

dl = (max(wzrost)-min(wzrost))/(k-1)

podzial2 = dl*(0:k) - dl/2 + min(wzrost)

szereg2 = table(cut(wzrost, podzial2))

szereg2

hist(wzrost, podzial2, col = "lightblue")

pstwa = pnorm(podzial2[2:(k+1)], mi, sigma) -pnorm(podzial2[1:k], mi, sigma)

pstwa

(t = chisq.test(szereg2, p = pstwa))

#ups, zle!

sum(pstwa)

#faktycznie, poprawimy tak

podzial2b = podzial2

podzial2b[1] = -Inf

podzial2b[length(podzial2b)] = Inf

pstwa = pnorm(podzial2b[2:(k+1)], mi, sigma) -pnorm(podzial2b[1:k], mi, sigma)

sum(pstwa)

(t = chisq.test(szereg2, p = pstwa))

#poprawiona p-wartość

1-pchisq(t$statistic, k-3)

#czegos jeszcze nie mamy - klasy po 5 elementow

szereg2 #uwaga to dziala dla k = 20,

podzial3 = podzial2b[c(1, 3:17, 21)]

(szereg3 = table(cut(wzrost, podzial3)))

k3 = length(szereg3)

pstwa3 = pnorm(podzial3[2:(k3+1)], mi, sigma) -pnorm(podzial3[1:k3], mi, sigma)

(t = chisq.test(szereg3, p = pstwa3))

#poprawiona p-wartość

1-pchisq(t$statistic, k3-3)

#to może prościej (i zarazem lepiej)? - test normalnosci Shapiro-Wilka

shapiro.test(wzrost)

#to jeszcze test niezależności chi^2

dane = read.csv2("Zakupy.csv")

attach(dane) #najprostszy test chi2, niezaleznosci

chisq.test(table(WYKSZTALCENIE, cut(WYDATEK, 7)))

#ale ostrzega (mało liczne klasy przy dużym wydatku)

#pierwsza poprawka

podzial = seq(min(WYDATEK), max(WYDATEK), length = 10)

table(cut(WYDATEK, podzial))

#no to tak:

podzial2 = podzial[c(1:7, 10)] #OK:

table(cut(WYDATEK, podzial2)) #to jeszcze to

wyk = factor(WYKSZTALCENIE, levels = c("P", "Z", "S", "W"))

table(wyk, cut(WYDATEK, podzial2))

chisq.test(table(wyk, cut(WYDATEK, podzial2)))

detach(dane)

Przetestować normalność dla danych z pliku wzrost.csv za pomocą testu chi-kwadrat, w oparciu o szereg rozdzielczy o 20 klasach o, w przybliżeniu, jednakowej liczności empirycznej.

wzrost=read.csv2("wzrost.csv")[,1]

wzrost

mi=mean(wzrost)

sigma=sd(wzrost)

w=seq(0,500,length.out=21)

podzial=quantile(wzrost, (w/500))

pstwa = pnorm(podzial[2:21], mi, sigma) - pnorm(podzial[1:20], mi, sigma)

# licznosci teoretyczne maja byc rowne dlatego to dajemy=09

szereg=table(cut(wzrost,podzial))

hist(wzrost,podzial)

barplot(szereg)

t =chisq.test(szereg, pstwa)

Przetestować normalność dla danych z pliku wzrost.csv za pomocą testu chi-kwadrat, w oparciu o szereg rozdzielczy o 15 klasach o jednakowej liczności teoretycznej (przykład z wykładu). Następnie wykonać ten sam test, ale przed utworzeniem szeregu rozdzielczego dokonać randomizacji danych, zaburzając je o wartość z rozkładu jednostajnego na przedziale [-½, ½].

wzrost =read.csv2("Wzrost.csv")[,1]

x=runif(500,-1/2,1/2)

wzrost=wzrost+x

mi=mean(wzrost)

sigma=sd(wzrost)

w=0:k1

k1=15

podzial =qnorm(w/k1, mi, sigma)

(szereg =table(cut(wzrost, podzial)))

k =length(szereg)

pstwa =pnorm(podzial[2:(k+1)], mi, sigma) -pnorm(podzial[1:k], mi,sigma)

(t =chisq.test(szereg, p =pstwa))

#poprawiona p-wartosc; to jest niepotrzebne bo licznosci teoretycznemamy rowne

1-pchisq(t$statistic, k-3)

Dla wszystkich (sześciu) par różnego wykształcenia (dane zakupy.csv) wykonać testy równości przeciętnego wydatku między grupami osób o danym wykształceniu.

zakupy=read.csv2("zakupy.csv")

attributes(zakupy$WYKSZTALCENIE)

p=zakupy$WYDATEK[(zakupy[,4]=="P")]

s=zakupy$WYDATEK[(zakupy[,4]=="S")]

w=zakupy$WYDATEK[(zakupy[,4]=="W")]

z=zakupy$WYDATEK[(zakupy[,4]=="Z")]

sd(p) sd(s) sd(w) sd(z)

t.test(p,s) t.test(p,w) t.test(p,z) t.test(s,w) t.test(s,z) t.test(w,z)

dane = read.csv2("Zakupy.csv")

attach(dane)

wyd = WYDATEK

plec = PLEC

wyksz = factor(WYKSZTALCENIE, levels = c("P", "Z", "S", "W"))

n =length(wyd)

detach(dane)

hist(wyd, 15, col = "green")

(poprzedni wykład)

by(wyd, wyksz, mean)

t.test(wyd[wyksz=="P"], wyd[wyksz=="S"])

t.test(wyd[wyksz=="P"], wyd[wyksz=="S"], var.equal = T)

#a gdyby trochę zmienić problem:

#przetestować hipotezę, że średnie wydatki we wszystkich grupach są takie same

#hipoteza alternatywna: istnieją dwie grupy w której się różnią

#lub inaczej: czy różnica którą zaobserwowaliśmy jest istotna statystycznie?

#trochę fałszywa (mediana, a nie średnie) ilustracja problemu:

boxplot(wyd~wyksz, col = c("blue", "grey", "green", "red"))

#no to poprawka - są średnie

points(1:4, by(wyd, wyksz, mean), pch=16, cex = 2)

#jedno rozwiązanie, niepozbawione wad, - test chi^2 niezależności

(sz_2w = table(cut2(wyd, g=5), wyksz))

chisq.test(sz_2w)

#zanim dojdziemy do celu - niby dygresja

#wariancja w ramach grup

by(wyd, wyksz, var)

#spróbujmy ją uśrednić - wariancja wewnątrzgrupowa

(w_wewn = sum(by(wyd, wyksz, var)*(by(wyd, wyksz, length)-1))/(n-1))

#to teraz coś drugiego - wariancja międzygrupowa

(w_miedz = sum(by(wyd, wyksz, length) * (by(wyd, wyksz, mean)-mean(wyd))^2)/(n-1))

#jak je dodamy

w_wewn + w_miedz

#to dostaniemy to:

var(wyd)

#intuicja mówi, że jeżeli podział na grupy daje małą wariancję

#międzygrupową, to średnie w ramach grup różnią się nieznacznie

#pierwszy na to wpadł sir Ronald Fisher

#bardziej formalnie, jeżeli zachodzi:

#1) niezależność zmiennych objaśniających (podziału na grupy) - tu nie ma problemu

#2) normalność rozkładu zmiennej objaśnianej w grupach - nie mamy

#3) jednorodność wariancji zmiennej objaśnianej w grupach - też raczej nie mamy

#(traktujmy więc poniższy przykład z ostrożnością)

#oraz prawdą jest równość średnich w grupach, to przy oznaczeniach:

#k to ilośc grup

k = 4

(MSTR = w_miedz * (n-1)/(k-1))

(MSE = w_wewn * (n-1)/(n-k))

#statystyka testowa, której wartość tu liczymy:

(stat_F = MSTR/MSE)

#ma rozkład Fishera-Snedecora o (k-1),(n-k) stopniach swobody

#zatem p-wartośc jest równa

1-pf(stat_F, k-1, n-k)

#czyli na poziomie 1-alfa = 0,95 nie mamy podstaw do odrzucenia hipotezy

#o równości średnich w grupach

#rzecz jest znana statystykom jako analiza wariancji (ANOVA)

#więc na szczęście w R nie trzeba liczyć tego na piechotę:

anova(lm(wyd~wyksz))

#to teraz kolejny przykład - sztuczny podział "płciowykształcenie"

c(wyksz)-1

(tmp_pw = c(wyksz)-1 + 4 * (c(plec)-1))

tmp_pw = factor(tmp_pw)

levels(tmp_pw) = c("KP", "KZ", "KS", "KW", "MP", "MZ", "MS", "MW")

boxplot(wyd~tmp_pw, col = heat.colors(8))

points(1:8, by(wyd, tmp_pw, mean), pch=16, cex = 2)

#anova

anova(lm(wyd~tmp_pw))

#to samo prościej

anova(lm(wyd~plec:wyksz))

#i oczywiście

boxplot(wyd~plec:wyksz, col = heat.colors(8))

points(1:8, by(wyd, plec:wyksz, mean), pch=16, cex = 2)

#ups trochę inaczej - boxplot i by dają inną kolejność

boxplot(wyd~wyksz:plec, col = heat.colors(8))

points(1:8, by(wyd, plec:wyksz, mean), pch=16, cex = 2)

#wracając do wyniku:

anova(lm(wyd~plec:wyksz))

#naturalne pytanie po odrzuceniu hipotezy - które grupy mają rózne średnie

#odpowiedź - testy post-hoc

#np. test t-Studenta, ale lepiej poprawić poziom istotności (test Bonferroniego)

# alfa' = 2 * alfa / k / (k-1) (bo testuję (k*(k-1)/2) par)

# jeszcze lepiej test HSD (Tukey'a)

TukeyHSD(aov(wyd~plec:wyksz))

plot(TukeyHSD(aov(wyd~plec:wyksz)))

#to może prostszy przykład - gdybyśmy testowali na poziomie 0,9

anova(lm(wyd~wyksz))

TukeyHSD(aov(wyd~wyksz), conf.level = 0.9)

plot(TukeyHSD(aov(wyd~wyksz), conf.level = 0.9))

#przykład wcale nie był prostszy - odrzuciliśmy równość wszystkich średnich

#ale nie umiemy wskazać średniej, która odstaje (wszystkie p > 0,1)

#to prawie na koniec - najprostszy przykład jednoczynnikowej anova:

anova(lm(wyd~plec))

#dla porównania test t-Studenta

t.test(wyd[plec=="K"], wyd[plec=="M"], var.equal = T)

#grupowanie moze odbywać się po więcej niż jednej zmiennej (czynniku)

#wieloczynnikowa analiza wariancji

#dochodzą tu jednak dodatkowe problemy

#łatwiej jest wykonać obliczenia, gdy ilości elementów w grupach są jednakowe

#pewnie można też (tzn. R może) analizować gdy ilości są różne,

#nie wiem jednak co oznaczają wyniki

#zobaczmy więc, ile el. ma najmniejsza grupa:

(z=min(by(wyd, plec:wyksz, length)))

#ograniczymy więc wszystkie grupy do z=21 elementów

#zakładamy, że są one ustawione losowo,

#lub na wszelki wypadek wybieramy losową permutację

perm = sample(1:n, n)

#nowa ramka danych

ndane = data.frame(wyd=wyd[perm], plec=plec[perm], wyksz=wyksz[perm])

#podział wg grup

podzial = split(ndane, ndane$plec:ndane$wyksz)

#ograniczam rozmiar do z=21 elementów

for (i in 1:8) {

podzial[[i]] = podzial[[i]][1:z,]

}

#teraz małe "czary", ale zadziała - połączenie po podziale

ndane = unsplit(podzial, rep(1:8, z))

#jak widać jest podobnie

summary(ndane)

summary(wyd)

#trochę się jednak zmieniło

anova(lm(ndane$wyd~ndane$plec))

anova(lm(ndane$wyd~ndane$wyksz))

anova(lm(ndane$wyd~ndane$plec:ndane$wyksz))

#teraz najważniejsze - przy wielu czynnikach - interakcje

#to, żeby dwa wykresy w jednym oknie, następnie zmiana marginesów

par(mfrow=c(2,1))

par(mar=c(4.2, 4, 0.8, 1.1))

with(ndane, interaction.plot(plec, wyksz, wyd, xlab=""))

#to samo, tylko w drugą stronę

with(ndane, interaction.plot(wyksz, plec, wyd, xlab=""))

#analiza wariancji bez uwzględnienia interakcji

with(ndane, anova(lm(wyd~plec+wyksz)))

#z uwzględnieniem interakcji

with(ndane, anova(lm(wyd~plec*wyksz)))

#to samo inaczej

with(ndane, anova(lm(wyd~plec+wyksz+plec:wyksz)))

with(ndane, TukeyHSD(aov(wyd~plec*wyksz)))

#co jeśli faktycznie brak normalności zmiennej objaśnianej w grupach

#można użyć testu Kruskala-Wallisa, który testuje równość median w grupach

#ale założenia są też silne - zbliżony kształt i skala rozkładów

#(tzn. dopuszczone tylko przesunięcie (?))

boxplot(wyd~wyksz, col = terrain.colors(4))

kruskal.test(wyd~wyksz)

#działa tylko wersja jednoczynnikowa, ale

#sztuczny podział "płciowykształcenie"

c(wyksz)-1

(tmp_pw = c(wyksz)-1 + 4 * (c(plec)-1))

tmp_pw = factor(tmp_pw)

levels(tmp_pw) = c("KP", "KZ", "KS", "KW", "MP", "MZ", "MS", "MW")

boxplot(wyd~tmp_pw, col = rainbow(8))

kruskal.test(wyd~tmp_pw)

Podzielić dane z pliku zakupy.csv na trzy mniej więcej jednakowo liczne grupy ze względu na wiek klienta (powiedzmy Y, M, O). Następnie zbadać, czy różnica średnich wydatków w tych grupach jest istotna statystycznie.

zakupy=read.csv2("zakupy.csv")

attach(zakupy)

wyd =WYDATEK

plec =PLEC

wiek=WIEK

wyksz =factor(WYKSZTALCENIE, levels =c("P", "Z", "S", "W"))

n =length(wyd)

detach(zakupy)

max(zakupy[,3])

min(zakupy[,3])

mi=mean(zakupy[,3])

sigma=sd(zakupy[,3])

w=seq(0,max(zakupy[,3]),length.out=4)

podzial=quantile(zakupy[,3], (w/59))

podzialb=podzial

podzialb[1]=19

podzialb

szereg=table(cut(zakupy[,3],podzialb))

szereg

table(wyd, cut(wiek, podzialb))

chisq.test(table(wyd, cut(wiek, podzialb)))

szereg2 =table(cut2(wiek, g=3))

table(cut2(wyd,g=10), cut2(wiek, g=3))

chisq.test(table(cut2(wyd,g=10), cut2(wiek, g=3)))

#niezaleznosc za slaba, trzeba skorzystac z czegos mocniejszego wiec z tego:

lm(wyd~cut2(wiek,g=3))

anova(lm(wyd~cut2(wiek,g=3)))

cat("Nie ma podstaw do odrzucenia hipotezy H0, czyli srednie nie SA istotne")

Wykonać zadanie podobne do poprzedniego, z tym, że podział na grupy uwzględniać ma także płeć klienta, dając w ten sposób 6 grup (powiedzmy YW, MW, OW, YM, MM, OM).

szereg3 =table(cut2(wiek, g=3),plec)

chisq.test(table(wyd, cut2(wiek, g=3),plec))

lm(wyd~cut2(wiek,g=3):plec)

anova(lm(wyd~cut2(wiek,g=3):plec))

cat("nie ma podstaw do odrzucenia hipotezy H0, czyli roznica srednich w tych grupach jest nie istotna statystycznie")

lm(cut2(wyd,g=9)~cut2(wiek,g=3):plec)

wieloczynnikowa analiza wariacji

library("Hmisc")

dane = read.csv2("anova2.csv")

dane

# stawiamy hipotezy zerowe:

#H_oi: Srednie czasy spadania klocków pokrytych danym materia³em

#H_1i: �rednie tarcie nie jest równe

#H_oii: Srednie dok³adno�ci pomiarów kazdego z uczniów sa równe.

#H_1ii: Srednie dokladno�ci pomiarów kazdego z uczniów nie s¹ równe

#H_oiii: Nie zachodzi interakcja miedzy cechami I i II

#H_1iii: Zachodzi interakcja miedzy cechami I i II lub próba nie jest losowa lub

# wyniki prób dla cech I i II zale¿¹ od innej cechy

#Z regu³y najpierw powinni�my sprawdziæ czy zachodzi interakcja miêdzy cechami

# I i II. Je�li nie to mo¿emy wykonywaæ dalsze testy.

#Je�li interakcja zachodzi (odrzucili�my hipotezê H_oiii) to

#badanie anova nie ma senu, a wyniki dla II i I cechy zale¿¹ np od innej cachy, lub próba nie by³a losowa.

#Wiêc przeprowadzamy najpierw test z interakcjami

with(dane, interaction.plot(material, uczen, czas))

# lub to:

with(dane, interaction.plot(uczen, material, czas))

anova(aov(czas~material*uczen,dane))

# Pr=0.08914-prawdopodobienstwo ¿e hipoteza zerowa jest fa³szywa jest niskie:

#nie mamy podstaw do odrzucenia hipotezy zerowej, ¿e nie ma interakcji miêdzy

#cechami I i II.

#Mamy podstawy odrzuciæ hipotezê H_oi ze �rednie czasy spadania

#klocków pokrytych danym materia³em s¹ równe

#Mamy podstawy odrzuciæ hipotezê H_oii,

# ¿e �rednie dok³adno�ci pomiarów ka¿dego z uczniów s¹ równe.

# Zbiór krytyczny dla interakcji I i II (czyli dla H_oiii) jest postaci

#WIII(alpha)=[x, +ininity)=[F_(p-1)*(r-1),N-rp)(alpha),+infinity)

# F_(p-1)(r-1),N-rp)(alpha) oznacza rozk³ad F Snedecora o (p-1)(r-1) stopniach swobody

#w liczniku i N-rp stopniach swobody w mianowaniku przy poziomie instotno�ci alpha

#gdy statystyka T nale¿y do zbioru krytycznego W to odrzucamy hipotezê zerow¹.

#-------------------------tylko dla przyk³adu --------------------------------------

#df(2.33,3,24 )=0.1 # 0.1 jest poziomem istotno�ci alpha=0.1,

#df(x,p-1=Df1,Df2=N-rp), gdzie x=F_(p-1,N-rp)(alpha)

# gdy TIII nale¿y do WIII to na badanym poziomie istotno�ci mamy podstawy odrzuciæ

# hipotezê zerow¹.

df(4, 6,24 ) #=0.008

#U nas WIII_(0.008)=[4, +inifinity) dla alpha=0.008 TIII nie nale¿y

# do WIII- nie mamy podsaw odrzucenia hipotezy H_oiii. Zatem nie ma interakcji

#miêdzy cechami I i II.

df(2.33, 3,24) # = 0.1 zatem WI_(0.1)=[2.33, +infty) oraz TI=54.435 nale¿y do WI

# mamy zatem podstawy odrzuciæ hipotezê zerow¹ ¿e �redni czas spadania klocka

# nie zale¿y od rodzaju materi³u (odrzucamy hipoteze zerow¹)na poziomie istotno�ci

#alpa=0.1

df(4.32,2,24) #=0.1 W_II(0.1)=[4.32, + infty), TII=23,661 nale¿y do zbioru krytycznego,

#mamy wiêc podstwy odrzuciæ hipotezê ¿e uczniowie dokonuja pomiarów z �rednio tak¹

#sam¹ dok³adno�ci¹ (odrzucamy H_oii) na poziomie istotno�ci (0.1)

#------------opis oznaczeñ----------------------------------------------------------

#Df to stopnie swobody

#Mean Sq oznaczaja wariancje MSI(wariancja kolumn) - materialu,

#MSII (wariancja wierszy)- ucznia,

# Mean Sq i Residuals= MSE - wariancja wewn¹trzgrupowa(b³êdu)

# uczen:material to MS(I,II)- czyli wariancja interakcji

# Fvalue to warto�ci poszczególnych testów TI=54.435=MSI/MSE,

#TII=23,661=MSII/MSE, TIII=2.113=MS(I,II)/MSE

# gdzie kolejne wariancje wyznaczamy przez podzielanie Sum Sq/Df

anova(aov(czas~material+uczen,dane)) # test anova bez interakcji

Dla danych z pliku wzrost_zal.csv zbadać ,jak kolor oczu i kolor wlosow wplywaja na sredni wzrost czlowieka(z badanej populacji). Narysować wykresy ramkowe (uwaga:przedstawiając mediane) oraz oba wykresy interakcji. Dopuscic rozna ilość obserwacji w grupach.

dane=read.csv2("wzrost_zal.csv")

attach(dane)

oczy=kolor.oczu

wlosy=kolor.włosów

wzrost=wzrost..cm.

detach(dane)

attributes(dane)

#z iteracja:

with(dane, anova(lm(wzrost~oczy*wlosy))) #albo

with(dane, anova(lm(wzrost~oczy*wlosy))) #to jest bez iteracji:

anova(lm(wzrost~wlosy*oczy))

boxplot(wzrost~oczy, col =terrain.colors(4))

boxplot(wzrost~wlosy, col =terrain.colors(4))

boxplot(wzrost~wlosy*oczy, col = terrain.colors(4))

with(dane,interaction.plot(oczy,wlosy,wzrost,xlab=""))

with(dane,interaction.plot(wlosy,oczy,wzrost,xlab=""))

Wykonac poprzednie zadanie ograniczajac (w sposob losowy) probke tak, aby we wszystkich grupach bylo tyle samo obserwacji.

#mieszamy nasze dane

z=min(by(wzrost,oczy:wlosy,length)) #nowa ramka danych

set.seed(1) #tak zeby miec takie same wyniki jak Pati

perm =sample(1:150, 150)

ndane =data.frame(wzrost=wzrost[perm], oczy=oczy[perm], wlosy=wlosy[perm])

podzial = split(ndane, ndane$oczy:ndane$wlosy)

for(i in 1:8) {

podzial[[i]]=podzial[[i]][1:z,]}

#CZyli laczymy wszystkie dane w jedna kolumne

ndane2=unsplit(podzial,rep(1:8,z))

anova(lm(ndane2$wzrost~ndane2$wlosy*ndane2$oczy))

Dla danych z pliku czas_wykonania.csv zbadac, jak zdobyte kwalifikacje oraz plec pracownika wplywaja na czas wykonania czynnosci.

CZAS=read.csv2("czas_wykonania.csv ")

with(CZAS, anova(lm(czas~kwalifikacje*plec)))

#Tutaj pmiedzy plcia a kwalifikacja czy zachodzi interakcja..

with(CZAS, anova(lm(czas~plec*kwalifikacje)))

#a tutaj dziele na szesc grup i widac czy mezczyzna o kwalifikacji 3 ma wplyw na wykonywanie czas..

with(CZAS, anova(lm(czas~plec:kwalifikacje)))

with(CZAS,interaction.plot(kwalifikacje,plec,czas,xlab=""))

with(CZAS,interaction.plot(plec,kwalifikacje,czas,xlab=""))

Dla danych zapisanych w pliku raty.csv:

Czy modele regresji uzyskane w zadaniu 1 różnią się w sposób istotny statystycznie (alfa = 0,05)?

Utworzyć liniowe modele regresji opisujące zależność wydatku od dochodów klienta z uwzględnieniem decyzji o zakupie ratalnym:

* przyjmując, że linie regresji, dla różnych decyzji o ratach, są równoległe

* nie przyjmując powyższego założenia

Dla obu modeli narysować wykresy linii regresji (oraz punktów danych)

Rozwiazanie dla podpunktu 1) w=aD+bR+c

dane=read.csv2("raty.csv")

attach(dane)

wydatek=Wydatek

dochody=Dochody

wiek=Wiek

raty=Raty

detach(dane)

#kropki

plot(dochody,wydatek,pch=c(1,1)[raty],col=c("green","blue")[raty])

reg1=lm(wydatek~dochody+raty)

ndane =data.frame(dochody = c(min(dochody), max(dochody)),raty=c("N","N"))

lines(ndane$dochody, predict(reg1, ndane), lwd = 2, col = "red")

ndane1 =data.frame(dochody =c(min(dochody), max(dochody)),raty=c("T","T"))

lines(ndane1$dochody, predict(reg1, ndane1), lwd = 2, col = "blue")

#podpunkt 2)

ndane=data.frame(dochody=c(min(dochody),max(dochody)))

plot(dochody,wydatek,pch=c(1,1)[raty],col=c("green","blue")[raty])

reg=lm(wydatek~dochody)

raty c(raty)

raty1=c(raty)-1

(dochR=dochody[raty1==1])

(dochBR=3Ddochody[raty1==0])

(wydatekR=wydatek[raty1==1])

(wydatekBR=3Dwydatek[raty1==0])

(reg1=3Dlm(wydatekR~dochR))

reg2=3Dlm(wydatekBR~dochBR)

ndaneR=data.frame(dochR=c(min(dochody),max(dochody)))

ndane2=data.frame(dochBR=c(min(dochody),max(dochody)))

#aby dorysowac linie do kropek nalezy uzyc funkcje lines

lines(ndaneR$dochR,predict(reg1,ndaneR),lwd=2,col="blue")

lines(ndane2$dochBR,predict(reg2,ndane2),lwd=2,col="green")

lines(ndane$dochody,predict(reg,ndane),lwd=2,col="black")

Utworzyć (i narysować wykres) model regresji liniowej opisujący zależność wydatku od (równocześnie) dochodów i wieku klienta

reg3=lm(wydatek~dochody+wiek)

reg3

summary(reg3)

(min(dochody))

(max(dochody))

(min(wiek))

(max(wiek))

(mv=seq(1000, 4200,by=450))

(ov=seq(21,60))

#funkcja regresji

f=function(x,y)

reg3$coef[1]+x*reg3$coef[2]+y*reg3$coef[3]

#macierz wartosci funkcji regresji:

dv=outer(mv,ov,f)

wykres=persp(x=mv, y=ov, z=dv, xlab="dochody", ylab="wiek", zlab="wydatek",ticktype="detailed",theta=36,phi=28)

points(trans3d(dochody, wiek, wydatek, wykres), col="red", lwd = 2, pch = 20)

#nanosimy dane empiryczne:

#oraz kreski empiryczne wyznaczajace polozenie punktow wzgledem =

plaszczyzny regresji:

segments(trans3d(dochody, wiek, wydatek, wykres)$x,

trans3d(dochody, wiek, wydatek, wykres)$y,

trans3d(dochody, wiek, predict(reg3), wykres)$x,

trans3d(dochody, wiek, predict(reg3), wykres)$y)

Utworzyć (i narysować wykres) model regresji opisujący decyzję o zakupie na raty w zależności od wieku klienta.

raty1=c(raty)-1

reg5=lm(raty1~wiek)

reg5

plot(wiek,raty1)

ndane=data.frame(wiek=c(min(wiek),max(wiek)))

ndane

lines(ndane$wiek, predict(reg5,ndane),lwd=3,col="green")

Utworzy liniowe modele regresji opisujace zaleznosc wydatku od dochodow klienta

z uwzglednieniem decyzji o zakupie ratalnym:

* przyjmujac, ze linie regresji, dla rnych decyzji o ratach, s rwnolege

* nie przyjmujc powyszego zaoenia

Dla obu modeli narysowa wykresy linii regresji (oraz punktw danych)

#Rozwizanie dla podpunktu 1) w=aD+bR+c

(dane=read.csv2("raty.csv"))

attach(dane)

wydatek=Wydatek

dochody=Dochody

wiek=Wiek

raty=Raty

detach(dane)

#punkty danych

plot(dochody,wydatek,pch=c(1,1)[raty],col=c("green","blue")[raty])

reg=lm(wydatek~dochody)

#Inny sposob nanoszenie linii regesji na wykres

#nanosimy krzywa regresji na wykres:

abline(reg,col="blue",lwd=2)

#nowa ramka danych, wystarczy opisac dwa punkty bo chcemy miec prosta i uwzgledniamy decyzje

(ndane = data.frame(dochody = c(min(dochody), max(dochody)),raty=c("N","N")))

#Aby narysowac wykres na punktach danych stosujemy polecenie "lines"

lines(ndane$dochody, predict(reg1, ndane), lwd = 2, col = "red")

#To samo postepowanie ale w przypadku decyzji na TAK

(ndane1 = data.frame(dochody = c(min(dochody), max(dochody)),raty=c("T","T")))

lines(ndane1$dochody, predict(reg1, ndane1), lwd = 2, col = "blue")

podpunkt 2)

(ndane2=data.frame(dochody=c(min(dochody),max(dochody))))

plot(dochody,wydatek,pch=c(1,1)[raty],col=c("green","blue")[raty])

(reg=lm(wydatek~dochody))

raty

c(raty)

(raty1=c(raty)-1)

#1==T, 0==N

(dochR=dochody[raty1==1])

(dochBR=dochody[raty1==0])

(wydatekR=wydatek[raty1==1])

(wydatekBR=wydatek[raty1==0])

(reg1=lm(wydatekR~dochR))

(reg2=lm(wydatekBR~dochBR))

(ndaneR=data.frame(dochR=c(min(dochody),max(dochody))))

(ndaneBR=data.frame(dochBR=c(min(dochody),max(dochody))))

#punkty danych

plot(dochody,wydatek,pch=c(1,1)[raty],col=c("green","blue")[raty])

lines(ndaneR$dochR,predict(reg1,ndaneR),lwd=2,col="blue")

lines(ndaneBR$dochBR,predict(reg2,ndaneBR),lwd=2,col="green")

#ta lini jaest niepotrzebna

lines(ndane$dochody,predict(reg,ndane),lwd=2,col="black")

Utworzy (i narysowa wykres) model regresji liniowej opisujcy zaleno wydatku

od (rwnoczenie) dochodw i wieku klienta

(reg3=lm(wydatek~dochody+wiek))

summary(reg3)

(min(dochody))

(max(dochody))

(min(wiek))

(max(wiek))

(mv=seq(1000, 4200,by=450))

(ov=seq(21,60))

#funkcja regresji

f=function(x,y)

reg3$coef[1]+x*reg3$coef[2]+y*reg3$coef[3]

#macierz wartosci funkcji regresji:

(dv=outer(mv,ov,f))

wykres=persp(x=mv, y=ov, z=dv, xlab="dochody", ylab="wiek", zlab="wydatek",ticktype="detailed",theta=36,phi=28)

points(trans3d(dochody, wiek, wydatek, wykres), col="red", lwd = 2, pch = 20)

#nanosimy dane empiryczne:

#oraz kreski empiryczne wyznaczajace polozenie punktow wzgledem plaszczyzny regresji:

segments(trans3d(dochody, wiek, wydatek, wykres)$x,

trans3d(dochody, wiek, wydatek, wykres)$y,

trans3d(dochody, wiek, predict(reg3), wykres)$x,

trans3d(dochody, wiek, predict(reg3), wykres)$y)

Utworzy (i narysowa wykres) model regresji opisujcy decyzj

o zakupie na raty w zalenoci od wieku klienta.

raty1=c(raty)-1

(reg5=lm(raty1~wiek))

plot(wiek,raty1)

(ndane3=data.frame(wiek=c(min(wiek),max(wiek))))

lines(ndane3$wiek, predict(reg5,ndane3),lwd=3,col="green")

Dla danych z pliku test_regresji.csv dopasowa krzyw obrazujc zaleno y do x

poprzez model regresji (parametrycznej) wielomianem stopnia 2.

Narysowa wykres.

(dane=read.csv2("test_regresji.csv"))

attributes(dane)

x=dane$x

y=dane$y

plot(x,y)

#nieliniowy model regresji tworzy sie za pomoca "nls"

reg=nls(y~a*x^2+b*x+d,start=list(a=0,b=0,d=0))

dane1=data.frame(x=seq(min(x),max(x),length=150))

lines(dane1$x,predict(reg,dane1),lwd=1,col="blue")

Dla powyszych danych dopasowa krzyw metod regresji lokalnie wielomianowej,

z uyciem domylnych opcji, a nastpnie zaobserwowa wpyw parametrw

sterujcych (span i degree) na sposb dopasowania. Narysowa odpowiednie wykresy.

#regresja lokalnie wielowymiarowa

(reg2=loess(y~x))

plot(x,y)

#span=1 bierze wszystkie dane span=0,75 bierze 0,75 danych- jest to stopein nieliniowosci

lines(dane1$x,predict(reg2,dane1),lwd=3,col="green")

lines(dane1$x,predict(loess(y~x,span=0.1, degree=1),dane),lwd=2,col="pink")

cat("Im mniejszy span tym linia przechodzi przez wiecej punktw,bardziej dopasowana do punktw")

lines(dane1$x,predict(loess(y~x,span=1, degree=1),dane1),lwd=2,col="grey")

lines(dane1$x,predict(loess(y~x,span=0.1, degree=0),dane1),lwd=2,col="purple")

lines(dane1$x,predict(loess(y~x,span=0.5, degree=1),dane1),lwd=2,col="orange")

#degree=0 malo uzywac

Dla tych samych danych dopasowa krzyw metod gadkich funkcji sklejanych

z uyciem domylnych opcji. Nastpnie zaobserwowa wpyw parametru spar

na sposb dopasowania. Narysowa odpowiednie wykresy

plot(x,y)

#gladkie funkcje sklejane, spar parametr wygladzania

reg3=smooth.spline(x,y,spar=0.1)

lines(predict(reg3,x=seq(min(x),max(x),length=4000)),lwd=2,col="red")

reg4=smooth.spline(x,y,spar=0.5)

lines(predict(reg4,x=seq(min(x),max(x),length=4000)),lwd=2,col="black")

reg5=smooth.spline(x,y,spar=1)

lines(predict(reg5,x=seq(min(x),max(x),length=4000)),lwd=2,col="green")

cat("Im wiekszy spar{on moze byc w rpzedziale [0,1] } tym funkcja jest bardziej gadka")

Dla danych z pliku bezrobocie.csv (czas pozostawania bez pracy w zalenoci

od wyksztacenia, dla osb zarejestrowanych w 2009r w PUP nr 1, d)

zbada zaleno midzy wyksztaceniem a czasem pozostawania bez pracy.

Narysowa map percepcji.

(dane2=read.csv2("Bezrobocie.csv",header=T))

names(dane2)

chisq.test(dane2)

#mamy zaleznosc

#zrobimy najpierw to rcznie

#macierz korespodencji(czstoi)

#bierzemy sum(dane2) poniewaz mamy gotowa tabelke,a na wykladzie musielismy wydobyc tabele

P=dane2/sum(dane2)

sum(P)

#sumowanie elementw w P powinno dac jeden

liczb_w=nrow(P)

liczb_k=ncol(P)

#masy wierszy i kolumn

(masa_w=rowSums(P))

(masa_k=colSums(P))

#Teoretyczne czstosci (w tescie chi-kwadrat)

teor=outer(masa_w,masa_k,"*")

teor

#standaryzowane reszty Pearsonowskie

(E=(P-teor)/teor^(1/2))

heatmap(as.matrix(dane2))

$dekompozycja maceirzy E

(S=svd(E))

#Wartosc $d jest uporzadkowana malejaco i ona oznacza znaczenie

#standaryzowane wsprzdna wierzy i kolumn

#dziele kolumny u i v przez pierwiastkiz mas

X=diag(1/sqrt(masa_w))%*%S$u

Y=diag(1/sqrt(masa_k))%*%S$v

$wykres

plot(rbind(X[1:5,1:5],Y[1:6,1:5]),col="white",xlab="",ylab="",main="Wykres")

#Specjalnie jest na bialo zeby potem dopisac nazwy i zeby te punkty niemazaywaly obrazu

#i tutaj jest problem zeby nazwy dopasowac..nie wiem czemu??

text(X[1:5,1:5],levels(colnames(dane2)),col="blue")

text(Y[1:6,1:5],levels(dane2),col="red")

#To samo ale przy uzyciu funkcji, najpierw trzeba pakiet zainstalowac ca

tabela=table(dane2[,1],dane2[,2])

tabela

plot(ca(dane2),mass=T)

rownames(dane2)

summary(ca(dane2))

Dla danych z pliku centroid.csv dokonać grupowania metodą k-średnich (dla k = 10, czyli na 10 klas)  z użyciem metryki euklidesowej na standaryzowanych danych. Następnie narysować wykres oryginalnych danych z podziałem na klasy oraz uzyskanymi środkami klas

#metoda k-średnich

metraz=read.csv2("Centroid.csv")

attach(metraz)

#standaryzowanie danych

metraz2 = data.frame(osoby = metraz$OSOBY/sd(metraz$OSOBY), metraz = metraz$METRAZ/sd(metraz$METRAZ))

metraz2

podzial=kmeans(metraz2,10)

plot(metraz,col=podzial$cluster)

###points(METRAZ[podzial$cluster], pch=19,col="black")

sr1=podzial$center[,1]

sr2=podzial$center[,2]

nowe_sr1=sr1*sd(OSOBY)

nowe_sr2=sr2*sd(METRAZ)

nowecen=data.frame(nowe_sr1,nowe_sr2)

points(nowecen,pch=15,col="yellow")

Podzielić dane z pliku dane_centr.csv (pierwszych 500 wartości) na 6 klas metodą k-średnich (metryka euklidesowa po standaryzacji) oraz metodą grupowania wokół centroidów (metryka taksówkowa po standaryzacji).
Porównać uzyskane wyniki.

danece=read.csv2("danecen.csv")

attach(danece)

names(danece)

kol1=danece[1:500,1]

kol2=danece[1:500,2]

dane2=data.frame(kol1/sd(kol1),kol2/sd(kol2))

mks<-kmeans(dane2,6)

library(cluster)

#metoda grupowania wokol centroidow

p<-pam(dane2,6,metric="manhattan")

plot(dane2,pch=mks$cluster,col=p$cluster)

points(mks$center,pch=17,col="black")

points(p$medoids, pch=19, col="yellow"cex=2)

Podzielić dane z pliku dane_centr.csv na 6 klas metodą grupowania woków centroidów (metryka euklidesowa po standaryzacji). Warto użyć wersji metody odpowiedniej dla dużej ilości danych.

attach(danece)

duze1=data.frame(x/sd(x),y/sd(y))

duze1

p2<-clara(duze1,6,metric="euclidean")

plot(duze1,col=p2$cluster)

points(p2$medoids,pch=19,col="black",cex=2)

e

Zad e

dane=read.csv2("dane_usmiech2.csv")

attach(dane)

dane2=data.frame(x=dane$x[klasa==1])

dane2

plot(dane2)

attach(dane2)

k=20

mi=mean(dane2)

sigma=sd(dane2)

ile=length(dane2)

w=0:k

podzial=qnorm(w/k,mi,sigma)

podzial

podzial[1]=min(dane2)-0.1

podzial[length(podzial)]=max(dane2)+0.1

szereg=table(cut(dane2[,],podzial))

szereg

(t=chisq.test(szereg))

t

#nie ma podstaw do odrzucenia

1-pchisq(t$statistic, k-3)

(hist(dane2[,], podzial, col = "yellow"))

attr(,"class")

plot(density(dane2[,]))

detach(dane2)

detach(dane)

Zad k

dane =read.csv2("regresja2.csv")

dane

attach(dane)

#regresja paramteryczna

reg = nls(y ~ a*(sin(b*x))+c , start=list(a=1,b=1,c=0))

summary(reg)

#regresja nie parametryczna

reg2=loess(y~x)

summary(reg2)

#wykresy

plot(x,y)

ndane = data.frame(x = seq(min(x), max(x), length = 300))

#parametryczny

lines(ndane$x, predict(reg, ndane), lwd = 2, col = "blue")

#nie parametryczny

lines(ndane$x, predict(reg2, ndane), lwd = 2, col = "red")

#obliczamy R_kwadrat dla regresji parametrycznej

a = 2.008283

b=1.001095

c = 0.014549

x=dane$x

y = a*(sin(b*x))+c

y_t=y

y=dane$y

roznica1= (y-mean(y))^2

roznica2 = (y_t-mean(y))^2

R_kwadrat = sum(roznica2)/sum(roznica1)

R_kwadrat

Zadanie 6.1

dane=read.csv2("zakupy.csv")

# wykonam test chi^2 niezaleznosci

attach(dane)

chisq.test(table(WYDATEK,(cut(WIEK,3)))) #jest ok, ale klasy mało liczne

podzial=seq(min(WYDATEK),max(WYDATEK), length=5)

table(cut(WYDATEK,podzial))

wiek=cut(WIEK,3)

levels(wiek)=c("Y", "M", "O")

table(wiek,cut(WYDATEK,podzial))

chisq.test(table(wiek,cut(WYDATEK,podzial))) #jeszcze można by pokombinować z gabarytami klas, ale już jest znośnie.

#2. Dopasować możliwie dobry model regresji (parametrycznej)

#do zależności y od x w 4-tej klasie danych dane_usmiech2.csv

dane=read.csv2("dane_usmiech2.csv")

dane

attach(dane)

x=dane$x[klasa=="4"]

y=dane$y[klasa=="4"]

x

y

plot(x,y)

reg = nls(y ~ a*x^2 + b*x + c , start=list(a=1,b=1,c=0))

summary(reg)

ndane = data.frame(x = seq(min(x), max(x), length = 300))

lines(ndane$x, predict(reg, ndane), lwd = 2, col = "blue")

dobra czyli trza zapamietac ze jak w zadaniu jest napisane ze "jednakowa licznosć teoretyczna"
to podzial liczymy w ten sposób:
k = 15
w = 0:k
mi = mean(wzrost)
sigma = sd(wzrost)
podzial = qnorm(w/k,mi, sigma)






a jak jest "jednakowa liczność empiryczna" to :
mi=mean(wzrost)
sigma=sd(wzrost)
w=seq(0,500,length.out=21)
w
podzial=quantile(wzrost, (w/500))