#-----------------------------------------------
# exemple de résolution matricielle du problème
# de la régression linéaire "à la main" 
# au lieu d'utiliser lm()
# P. Coquillard 2019
#-----------------------------------------------

### initialisation du générateur aléatoire 
  set.seed(241)

### On simule des données
  nobs <- 250  # nombre observations
  b0 <- 4      # le futur intercept
  b1 <- 2      # le futur coef directeur

# maintenant on génère 250 données fictives
# en ajoutant un peu de bruit avec rnorm()
# qui par défaut tire dans N(0,1)

  x <- rnorm(nobs)
  y <- b0 + b1*x + rnorm(nobs, 0, 0.5)
  df <- data.frame(x, y)

# voyons le graphique des données (scatter plot)
  plot(x,y)

### la solution analytique

# d'abord construire la matrice X
  X <- model.matrix(y ~ x, data = df)

# resoudre l'équation matricielle 
# b = (X'X)^-1 * X'Y comme vu en cours
  beta <- solve(t(X) %*% X) %*% t(X) %*% y
  beta
  abline(beta)

# les betas obtenus par régression sont très proches des
# betas utilisés pour fabriquer les données fictives
# NB : dans le cas d'une régression multiple, la démarche
# aurait été exactement la même

#  vérifions  
  fit <- lm(y~x)
  summary(fit)

# CQFD !  ;-)

###======================EOF===========================###