#### Lineární modely - BLUE #####

# databáze: 4 dojnice: 4000, 5000, 6000, 7000 kg mléka
# 1-2 dojnice ve stáde A, 3-4 dojnice ve stáde B
# databáze
#Jedinec c. Užitkovost (y) Stádo
# 1         4000            A
# 2         5000            A
# 3         6000            B
# 4         7000            B

############################
# príklad 1 - spocítat prumer pomocí lineárního modelu

# lineární model: y = mí + e

# maticový zápis: y = Xb + e

# soustava normálních rovnic: [X´X][b]=[X´y]

# matice X
X <- matrix(c(1,1,1,1),4,1)
X
# vektor y - užitkovosti
y <- matrix(c(4000,5000,6000,7000),4,1)
y
# X´X

XX <- t(X)%*%X
XX

Xy <- t(X)%*%y
Xy

b <- solve(XX) %*% Xy
b

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

# príklad 2 - zjistit prumery obou stád

# lineární model: y = S + e

# maticový zápis: y = Xb + e

# soustava normálních rovnic: [X´X][b]=[X´y]

X <- matrix(c(1,1,0,0,
              0,0,1,1),4,2)
X

XX <- t(X)%*%X
XX
Xy <- t(X)%*%y
Xy

b <- solve(XX) %*% Xy
b

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

# príklad 3 - soucasný odhad celkového prumeru a prumeru stád

# lienární model: y = mí + S + e
# v maticích: y = X1a + X2b + e
# soustava normálních rovnic: 
#  |X1´X1   X1´X2|    |a|    = |X1´y|
#  |X2´X1   X2´X2|    |b|      |X2´y|

# X1 - prumer celkový
X1 <- matrix(c(1,1,1,1),4,1)
X1

# X2 - prumer stád
X2 <- matrix(c(1,1,0,0,
               0,0,1,1),4,2)
X2

X1X1 <- t(X1)%*%X1
X1X1

X1X2 <- t(X1)%*%X2
X1X2

X2X1 <- t(X2)%*%X1
X2X1

X2X2 <- t(X2)%*%X2
X2X2

LS1 <- cbind(X1X1,X1X2) 
  LS1
LS2 <- cbind(X2X1,X2X2)
  LS2

LS <- rbind(LS1,LS2)
LS

X1y <- t(X1)%*%y
X1y

X2y <- t(X2)%*%y
X2y

PS <- rbind(X1y,X2y)
PS

det(LS)
  # determinant je roven nule -> matice LS je singulární a nejde ji invertovat  
 
# dodání podmínek rešitelnosti: prumer = 0, že stáda budou vyjádrena jako odchylky od prumeru

#     [,1] [,2] [,3]
# [1,]  4    2    2   |0|
# [2,]  2    2    0   |1|
# [3,]  2    0    2   |1|
#      |0    1    1|  |0|

a <- matrix(c(0,1,1),1,3)
a
b <- matrix(c(0,1,1,0),4,1)
b
LSa <- rbind(LS,a)
  LSa
LSab <- cbind(LSa,b)
  LSab
  
det(LSab)
  
  # matici PS musíme pridat jeden rádek s nulou

PS0 <- rbind(PS,0)
PS0

ab <- solve(LSab) %*% PS0
ab
