# 4. HARJOITUKSET
#################
# 4.1
# (a)
Y<-matrix(c(2,8,6,8,12,9,9,10),ncol=2)
S<-cov(Y)
y.ka<-cbind(apply(Y,2,mean))
m0<-cbind(c(7,11))
n<-nrow(Y); p<-ncol(Y)
T2<-n*t(y.ka-m0)%*%solve(S)%*%(y.ka-m0)
T2
# 13.63636

# (b)
# Hotellingin T2-jakaumaa vap. astein 2 ja 3

# (c)
F<-((n-p)/(p*(n-1)))*T2
F
# 4.545455
1-pf(F, p, n-p)  # Testisuureen arvoon liittyvä tn
          
# [1,] 0.1803279
# H0 jää voimaan

# 4.3
# (a)
Y<-read.table("C:\\Kurssit\\Mmm\\mmm03\\Datat\\lapset.txt", header = T, skip=5)
S1<-cov(Y[1:6,])
S2<-cov(Y[7:15,])
n1<-6; n2<-9
S<-((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)

S
          pituus    rinta   kvarsi
pituus 27.230769 6.561538 2.846154
rinta   6.561538 2.432308 1.400000
kvarsi  2.846154 1.400000 1.846154

# (b)
# H0: m1=m2
y1<-cbind(apply(Y[1:6,],2,mean))
y2<-cbind(apply(Y[7:15,],2,mean))
n<-n1+n2
p<-ncol(Y)
T2<-(n1*n2/n)*t(y1-y2)%*%solve(S)%*%(y1-y2)
T2
F<-((n-p-1)/(p*(n-2)))*T2
F
# 1.498179
 
1-pf(F, p,n1+n2-p-1)  # Testisuureen arvoon liittyvä tn
                
#  [1,]  0.2692616
          
# H0 jää voimaan
# Poikien ja tyttöjen keskiarvot eivät eroa merkitsevästi toisistaan
# (näin pienessä otoksessa)

# 4.4
# (b)
C<-matrix(c(1,0,0,0,0),nrow=4,ncol=5)
C1<-matrix(c(0,0,0,0,1),nrow=4,ncol=5)
C<-C-C1
C
     [,1] [,2] [,3] [,4] [,5]
[1,]    1   -1    0    0    0
[2,]    0    1   -1    0    0
[3,]    0    0    1   -1    0
[4,]    0    0    0    1   -1

ym<-cbind(c(20.8, 24.4, 22.6, 19.2, 14.0))
yn<-cbind(c(24.8,21.8,24.6,23.2, 20.4))
s<-c(72.7,33.0,71.7,18.7,22.3,0,21.3,21.3,12.7,11.9,0,0,41.3,16.4,9.8,0,0,0,11.5,10.2,0,0,0,0,21.7)
S<-matrix(s,nrow=5,ncol=5)
S<-S+t(S)-diag(diag(S))
S
d<-ym-yn
n1<-5; n2<-5; n<-n1+n2; p<-4
T2<-(n1*n2/n)*t(C%*%d)%*%solve(C%*%S%*%t(C))%*%C%*%d
T2
# 16.95028

F<-((n-p-1)/(p*(n-2)))*T2
F
# 2.648481
 
1-pf(F, p,n1+n2-p-1)  # Testisuureen arvoon liittyvä tn
                
#  0.1570882
          
#  Koska testisuureen arvoon liittyvä tn on suuri, niin H0 jää voimaan.

# 4.5
# (a) 
yv<-cbind(c(46.1,57.3,50.4))
S<-matrix(c(101.3,63,71,63,80.2,55.6,71,55.6,97.4),nrow=3)
C<-matrix(c(1,0,-1,1,0,-1),nrow=2,ncol=3)
 C
#     [,1] [,2] [,3]
#[1,]    1   -1    0
#[2,]    0    1   -1
n<-40
T2<-n*t(C%*%yv)%*%solve(C%*%S%*%t(C))%*%C%*%yv
T2
#[1,] 90.49458
p<-3
f<-((n-p+1)/((p-1)*(n-1)))*T2
F<-((n-p+1)/((p-1)*(n-1)))*T2
F
#[1,] 44.0871
1-pf(F, p-1,n-p+1)   # Testisuureen arvoon liittyvä tn

# 1.251548e-10
# H0 voidaan hylätä, riskitaso on lähellä nollaa
 
# (b) Parittaiset testit
# 5%:n riskitaso
n<-40
c1<-cbind(c(1,-1,0)); c2<-cbind(c(1,0,-1)); c3<-cbind(c(0,1,-1))
s1<-t(c1)%*%S%*%c1; s2<-t(c2)%*%S%*%c2; s3<-t(c3)%*%S%*%c3;
t<-qt(0.975,n-1)
t(c1)%*%yv-t*sqrt(s1/n); t(c1)%*%yv+t*sqrt(s1/n)
# t=2.382572      
# -13.58257  -8.817428

t(c2)%*%yv-t*sqrt(s2/n); t(c2)%*%yv+t*sqrt(s2/n) 
# -6.708192 -1.891808
# t=2.408192

t(c3)%*%yv-t*sqrt(s3/n); t(c3)%*%yv+t*sqrt(s3/n)
# 4.293945   9.506055
# t=2.606055        
# Yksittäiset hypoteesit hylätään tavallisilla t-testeillä testaten
# Todelliset 5%:n lv:t ovat leveämpiä. Kun testataan monta testiä peräkkäin samasta
# aineistosta, riskitaso ei ole enää sama kuin yksittäisessä testissä.

 
# 4.6
# (a)
Y<-read.table("C:\\Kurssit\\Mmm\\mmm03\\Datat\\pojat.txt", header = T, skip=7)
Y<-read.table("C:\\Kurssit\\Mmm\\mmm03\\Datat\\pojat.txt", header = T, skip=7)
Y<-transform(Y,x1=ppit1-ppit2,x2=plev1-plev2)
apply(Y[,5:6],2,mean)
#  x1   x2 
# 4.08 2.16
xv<-apply(Y[,5:6],2,mean)
S<-cov(Y[,5:6])
n<-25
T2<-n*t(xv)%*%solve(S)%*%xv
T2
# 8.541245     # T2(2,25)
p<-2; f<-n-1
F<-((f-p+1)/(f*p))*T2
1-pf(F, p,f-p+1)  # Testisuureen arvoon liittyvä tn
0.03016035
 
# Pään mittojen ero on merkitsevä 5%:n tasolla, mutta ei 1%:n tasolla

# 4.7
# (a)
plot(x1,x2)
points(xv[1],xv[2],col = "red")
points(0,0,col = "red")

# (b)
a<-solve(S)%*%xv
a<-solve(S)%*%xv
a
#           [,1]
#[1,] 0.05251437
#[2,] 0.05897740

#t(a)^2=T2, joten testaus palautuu 4.6a:n tilanteeseen