2.1 Opérateurs mathématiques

5 / 2

## [1] 2.5

5 %/% 2

## [1] 2

5 %% 2

## [1] 1

x <- c(2, 5, 7, 3)
y <- c(3, 5, 6, 4)

x > y

## [1] FALSE FALSE  TRUE FALSE

x != 5

## [1]  TRUE FALSE  TRUE  TRUE

FALSE < TRUE 

## [1] TRUE

caracteres_speciaux <- 
  c("!", "\"", "#", "$", "%", "&", "'", "(", ")", "*", "+", ",", "-", ".", "/", ":", ";", 
    "<", "=",">", "?", "@", "[", "\\", "]", "^", "_", "{", "|", "}", "~")
lettres_accentuees <- c("à", "â", "é", "è", "ê", "ë", "ï", "î", "ô", "ù", "ü", "û", "ç")
catacteres_ordonnes <- sort(c(caracteres_speciaux, 0:9, letters, LETTERS, 
                              lettres_accentuees, toupper(lettres_accentuees))) 
paste(catacteres_ordonnes, collapse = "")

"'-!\"#$%&()*,./:;?@[\\]^_{|}~+<=>0123456789aAàÀâÂbBcCçÇdDeEéÉèÈêÊëËfFgGhHiIîÎïÏjJkKlLmM
nNoOôÔpPqQrRsStTuUùÙûÛüÜvVwWxXyYzZ"

"arborescence" < "arbre" 

## [1] TRUE

"a" < "aa"

## [1] TRUE

2 < 10

## [1] TRUE

"2" < "10"

## [1] FALSE

!c(TRUE, FALSE)

## [1] FALSE  TRUE

p <- rep(c(FALSE, TRUE), each = 2)
q <- rep(c(FALSE, TRUE), times = 2)
cbind(p, q, "p et q" = p & q, "p ou q" = p | q, "p xor q" = xor(p, q))

##          p     q p et q p ou q p xor q
## [1,] FALSE FALSE  FALSE  FALSE   FALSE
## [2,] FALSE  TRUE  FALSE   TRUE    TRUE
## [3,]  TRUE FALSE  FALSE   TRUE    TRUE
## [4,]  TRUE  TRUE   TRUE   TRUE   FALSE

2 + 3 * 4

## [1] 14

(2 + 3) * 4

## [1] 20

2.2 Fonctions mathématiques opérant de façon vectorielle

R offre aussi plusieurs fonctions de calculs mathématiques, travaillant de façon vectorielle, dont les suivantes :

• racine carrée et fonctions relatives au signe : sqrt, abs, sign;
• exponentielles et logarithmes : exp, log (= logarithme naturel), log10, log2;
• fonctions trigonométriques : sin, cos, tan, acos, asin, atan, atan2;
• fonctions d’arrondissement : ceiling, floor, round, trunc, signif;
• fonctions reliées aux fonctions mathématiques bêta et gamma: beta, gamma, factorial, choose, etc.

Ces fonctions font un calcul distinct pour tous les éléments de l’objet fourni en entrée et retournent un résultat de même dimension que l’objet en entrée. Voici quelques exemples.

# Vecteur de données numériques pour les exemples
x <- seq(from = -1.25, to = 1.5, by = 0.25)
x

##  [1] -1.25 -1.00 -0.75 -0.50 -0.25  0.00  0.25  0.50  0.75  1.00  1.25  1.50

# Arrondissement régulier au dixième près
round(x, digits = 1)

##  [1] -1.2 -1.0 -0.8 -0.5 -0.2  0.0  0.2  0.5  0.8  1.0  1.2  1.5

# Arrondissement à l'entier supérieur
ceiling(x)

##  [1] -1 -1  0  0  0  0  1  1  1  1  2  2

# Arrondissement à la partie entière
trunc(x)

##  [1] -1 -1  0  0  0  0  0  0  0  1  1  1

Ces fonctions arrivent aussi à effectuer des calculs par élément dans un objet atomique de dimension supérieure à un ou dans un data frame.

# Matrice de données numériques pour les exemples
x_mat <- matrix(x, nrow = 2)
x_mat

##       [,1]  [,2]  [,3] [,4] [,5] [,6]
## [1,] -1.25 -0.75 -0.25 0.25 0.75 1.25
## [2,] -1.00 -0.50  0.00 0.50 1.00 1.50

# Extraction du signe
sign(x_mat)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]   -1   -1   -1    1    1    1
## [2,]   -1   -1    0    1    1    1

2.3 Fonctions mathématiques combinant des éléments

Certaines fonctions mathématiques en R effectuent des calculs faisant intervenir plus d’un élément de l’objet donné en entrée, plutôt que d’effectuer un calcul distinct pour chacun des éléments. C’est le cas des fonctions suivantes :

• somme ou produit de tous les éléments (retourne une seule valeur) : sum, prod;
• somme ou produit cumulatif des éléments (retourne un vecteur de même longueur que le vecteur en entrée) : cummsum, cumprod;
• différences entre des éléments : diff.

Voici quelques exemples.

# Matrice de données numériques pour les exemples
mat <- matrix(c(2, 5, 3, 4, 6, 5, 4, 3, 1, 2, 9, 8), nrow = 3, ncol = 4)
mat

##      [,1] [,2] [,3] [,4]
## [1,]    2    4    4    2
## [2,]    5    6    3    9
## [3,]    3    5    1    8

# Produit de tous les éléments
prod(mat)

## [1] 6220800

# Somme cumulative des éléments (ici 2, 2+5, 2+5+3, 2+5+3+4, ...)
cumsum(mat)

##  [1]  2  7 10 14 20 25 29 32 33 35 44 52

diff(mat)

##      [,1] [,2] [,3] [,4]
## [1,]    3    2   -1    7
## [2,]   -2   -1   -2   -1

mat[-1, ] - mat[-nrow(mat), ]

##      [,1] [,2] [,3] [,4]
## [1,]    3    2   -1    7
## [2,]   -2   -1   -2   -1

diff(c(2, 5, 3, 4))

## [1]  3 -2  1

diff(c(2, 5, 3, 4), lag = 2)

## [1]  1 -1

# soustractions effectuées : 3-2 et 4-5

2.4 Opérations sur des ensembles

Les fonctions R d’opérations sur des ensembles sont les suivantes :

• union : union,
• intersect : intersection,
• setdiff : différence,
• setequal : test d’égalité,
• is.element : test d’appartenance.

Les quatre premières fonctions permettent de comparer de différentes façons deux ensembles. Pour présenter des exemples d’utilisation des fonctions, créons donc d’abord deux ensembles fictifs. Ceux-ci sont stockés sous forme de vecteurs puisque le langage R de base n’offre pas de structure de données spécifique aux ensembles (contrairement au langage Python par exemple qui propose des sets).

A <- c("m", "s", "e", "f", "a")
B <- c("m", "e", "h", "i")

Rappelons que, par définition, un ensemble mathématique est une collection non ordonnée d’éléments uniques et, en conséquence, ne contient aucun élément dupliqué.

union(A, B)

## [1] "m" "s" "e" "f" "a" "h" "i"

c(A, B)

## [1] "m" "s" "e" "f" "a" "m" "e" "h" "i"

unique(c(A, B))

## [1] "m" "s" "e" "f" "a" "h" "i"

intersect(A, B)

## [1] "m" "e"

setdiff(A, B)

## [1] "s" "f" "a"

setdiff(B, A)

## [1] "h" "i"

setequal(A, B)

## [1] FALSE

setequal(c(5, 3, 2), c(2, 3, 2, 5, 5))

## [1] TRUE

is.element(el = c("d", "e"), set = A)

## [1] FALSE  TRUE

2.5 Calcul de distances

iris_ech_num <- iris[1:5, c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")]
iris_ech_num

##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1          5.1         3.5          1.4         0.2
## 2          4.9         3.0          1.4         0.2
## 3          4.7         3.2          1.3         0.2
## 4          4.6         3.1          1.5         0.2
## 5          5.0         3.6          1.4         0.2

dist(x = iris_ech_num, method = "euclidean", diag = TRUE)

##           1         2         3         4         5
## 1 0.0000000                                        
## 2 0.5385165 0.0000000                              
## 3 0.5099020 0.3000000 0.0000000                    
## 4 0.6480741 0.3316625 0.2449490 0.0000000          
## 5 0.1414214 0.6082763 0.5099020 0.6480741 0.0000000

dist(x = iris_ech_num, method = "minkowski", p = 2, diag = TRUE)

##           1         2         3         4         5
## 1 0.0000000                                        
## 2 0.5385165 0.0000000                              
## 3 0.5099020 0.3000000 0.0000000                    
## 4 0.6480741 0.3316625 0.2449490 0.0000000          
## 5 0.1414214 0.6082763 0.5099020 0.6480741 0.0000000

adist(x = "Allo", y = "Hello")

##      [,1]
## [1,]    2

2.6 Constantes mathématiques

En R, le nombre π est représenté par la constante pi.

pi

## [1] 3.141593

Inf est la constante R pour l’infini ∞.

-5 / 0

## [1] -Inf

NaN est une constante signifiant Not A Number. Cette constante est retournée par R lorsqu’un utilisateur lui demande d’effectuer une opération mathématique impossible, par exemple :

log(-1)

## Warning in log(-1): NaNs produced

## [1] NaN

Rappel : Attention à ne pas confondre la constante NaN avec la constante NA qui signifie plutôt Not Available et qui sert à représenter les données manquantes.

3.1 Conditions logiques vectorielles

Les deux premières utilités potentielles des conditions logiques énumérées ci-dessus requièrent la création d’un vecteur de valeurs logiques de la même longueur que l’objet R sur lequel la condition est testée. Nous avons vu au début de cette fiche des outils pour écrire de telles conditions logiques :

• les opérateurs de comparaison : ==, !=, >, >=, < et <=;
• les opérateurs et fonctions logiques vectoriels : ! (négation), & (et), | (ou) et xor (ou exclusif).

Voici des exemples d’écriture de conditions logiques utilisant le vecteur suivant, que nous avons déjà manipulé dans des notes précédentes.

de <- c(2, 3, 4, 1, 2, 3, 5, 6, 5, 4)

Supposons que nous voulions connaître le nombre d’éléments dans ce vecteur numérique dont la valeur est supérieure à 3. La condition logique suivante nous permet d’identifier ces valeurs.

condition <- de > 3
condition

##  [1] FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE

Compter le nombre de valeurs supérieures à 3 dans de revient à compter le nombre de TRUE dans le vecteur précédent. Ce calcul se réalise facilement avec la fonction sum comme suit.

sum(condition)

## [1] 5

Même si une somme est une opération mathématique sur des valeurs numériques, la commande précédente ne retourne par d’erreur, car R réalise d’abord une conversion implicite de type de données pour transformer les valeurs logiques en nombres (TRUE devient 1 et FALSE devient 0), puis effectue la somme.

de[condition]

## [1] 4 5 6 5 4

which(condition)

## [1]  3  7  8  9 10

de >= 3 & de <= 5

##  [1] FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE

de == 1 | de == 4 | de == 6

##  [1] FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE

!(de == 1 | de == 4 | de == 6)

##  [1]  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE

de != 1 & de != 4 & de != 6

##  [1]  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE

de %in% c(1, 4, 6)

##  [1] FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE

! de %in% c(1, 4, 6)

##  [1]  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE

c(1, 2, NA, 4, 5) == NA

## [1] NA NA NA NA NA

is.na(c(1, 2, NA, 4, 5))

## [1] FALSE FALSE  TRUE FALSE FALSE

is.finite(c(NA, NaN, Inf, -Inf, 4.5, 3))

## [1] FALSE FALSE FALSE FALSE  TRUE  TRUE

3.2 Conditions logiques de longueur 1

Lors de l’écriture d’une condition logique, il faut parfois s’assurer de retourner un vecteur logique de longueur 1. C’est le cas lors de l’écriture d’une condition logique dans une structure de contrôle conditionnelle if ... else (que nous verrons plus loin). La condition dans un if doit être obligatoirement de longueur 1.

de == 1 || de == 4 || de == 6

## [1] FALSE

is.logical(x) && length(x) == 1 && !is.na(x) && x   # isTRUE
is.logical(x) && length(x) == 1 && !is.na(x) && !x  # isFALSE

all(de %% 1 == 0)

## [1] TRUE

any(de < 0)

## [1] FALSE

is.numeric(de)

## [1] TRUE

5.1 Algèbre linéaire

Il existe plusieurs fonctions en R pour faire de l’algèbre linéaire.

• multiplication matricielle : %*%;
• transposition : t;
• inversion : solve (en fait solve résout A %*% x = B, mais par défaut B est la matrice identité);
• produit vectoriel (en anglais cross product) de matrices : crossprod;
• produit dyadique généralisé (en anglais outer product) : outer, %o%;
• produit de Kronecker généralisé : kronecker, %x%;
• matrices diagonales : diag;
• déterminant : det;
• valeurs et vecteur propres : eigen;
• décompositions : svd, qr, chol.

A <- matrix(1:6, nrow = 3, ncol = 2)
A

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6

B <- matrix(6:1, nrow = 3, ncol = 2)
B

##      [,1] [,2]
## [1,]    6    3
## [2,]    5    2
## [3,]    4    1

A*B

##      [,1] [,2]
## [1,]    6   12
## [2,]   10   10
## [3,]   12    6

A %*% B

## Error in A %*% B: non-conformable arguments

C <- matrix(c(5, 2, 3, 7), nrow = 2, ncol = 2)
C

##      [,1] [,2]
## [1,]    5    3
## [2,]    2    7

A %*% C

##      [,1] [,2]
## [1,]   13   31
## [2,]   20   41
## [3,]   27   51

solve(C)

##             [,1]       [,2]
## [1,]  0.24137931 -0.1034483
## [2,] -0.06896552  0.1724138

crossprod(A, B)

##      [,1] [,2]
## [1,]   28   10
## [2,]   73   28

# équivalent à
t(A) %*% B

##      [,1] [,2]
## [1,]   28   10
## [2,]   73   28

1:3 %o% 4:5

##      [,1] [,2]
## [1,]    4    5
## [2,]    8   10
## [3,]   12   15

1:3 %x% 4:5

## [1]  4  5  8 10 12 15

A <- matrix(12:1, nrow = 3, ncol = 4)
A

##      [,1] [,2] [,3] [,4]
## [1,]   12    9    6    3
## [2,]   11    8    5    2
## [3,]   10    7    4    1

B <- matrix(c(1,2), nrow = 2, ncol = 1)
B

##      [,1]
## [1,]    1
## [2,]    2

A %o% B

## , , 1, 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]   12    9    6    3
## [2,]   11    8    5    2
## [3,]   10    7    4    1
## 
## , , 2, 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]   24   18   12    6
## [2,]   22   16   10    4
## [3,]   20   14    8    2

A %x% B

##      [,1] [,2] [,3] [,4]
## [1,]   12    9    6    3
## [2,]   24   18   12    6
## [3,]   11    8    5    2
## [4,]   22   16   10    4
## [5,]   10    7    4    1
## [6,]   20   14    8    2

outer(1:3, 4:5, FUN = '+')

##      [,1] [,2]
## [1,]    5    6
## [2,]    6    7
## [3,]    7    8

kronecker(1:3, 4:5, FUN = '+')

## [1] 5 6 6 7 7 8

C

##      [,1] [,2]
## [1,]    5    3
## [2,]    2    7

diag(C)

## [1] 5 7

diag(1:3)

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

diag(3)

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

5.2 Calcul différentiel et intégral

(Ce qui est offert en R dans ce domaine n’est pas très performant ni facile d’utilisation.)

df <- deriv(expr = ~ log(x) + sin(x), namevec = "x")
df

## expression({
##     .value <- log(x) + sin(x)
##     .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
##     .grad[, "x"] <- 1/x + cos(x)
##     attr(.value, "gradient") <- .grad
##     .value
## })

x <- 2:5

eval(df)

## [1] 1.6024446 1.2397323 0.6294919 0.6505136
## attr(,"gradient")
##                x
## [1,]  0.08385316
## [2,] -0.65665916
## [3,] -0.40364362
## [4,]  0.48366219

log(x) + sin(x)

## [1] 1.6024446 1.2397323 0.6294919 0.6505136

1/x + cos(x)

## [1]  0.08385316 -0.65665916 -0.40364362  0.48366219

# Points en lesquels nous allons dériver
x <- as.double(-3:3)
# Valeur de la fonction en ces points
pnorm(x)

## [1] 0.001349898 0.022750132 0.158655254 0.500000000 0.841344746 0.977249868 0.998650102

# Calcul de la dérivée en ces points
numericDeriv(expr = quote(pnorm(x)), theta = "x")

## [1] 0.001349898 0.022750132 0.158655254 0.500000000 0.841344746 0.977249868 0.998650102
## attr(,"gradient")
##             [,1]       [,2]      [,3]      [,4]      [,5]       [,6]        [,7]
## [1,] 0.004431849 0.00000000 0.0000000 0.0000000 0.0000000 0.00000000 0.000000000
## [2,] 0.000000000 0.05399097 0.0000000 0.0000000 0.0000000 0.00000000 0.000000000
## [3,] 0.000000000 0.00000000 0.2419707 0.0000000 0.0000000 0.00000000 0.000000000
## [4,] 0.000000000 0.00000000 0.0000000 0.3989423 0.0000000 0.00000000 0.000000000
## [5,] 0.000000000 0.00000000 0.0000000 0.0000000 0.2419707 0.00000000 0.000000000
## [6,] 0.000000000 0.00000000 0.0000000 0.0000000 0.0000000 0.05399096 0.000000000
## [7,] 0.000000000 0.00000000 0.0000000 0.0000000 0.0000000 0.00000000 0.004431849

dnorm(x)

## [1] 0.004431848 0.053990967 0.241970725 0.398942280 0.241970725 0.053990967 0.004431848

delta <- .000001
(pnorm(x + delta) - pnorm(x - delta)) / (2 * delta)

## [1] 0.004431848 0.053990967 0.241970724 0.398942280 0.241970725 0.053990967 0.004431848

integrate(f = dnorm, lower = -Inf, upper = 1)

## 0.8413448 with absolute error < 1.5e-05

pnorm(1)

## [1] 0.8413447

5.3 Optimisation numérique

En mathématiques, l’optimisation consiste à trouver en quel(s) point(s) une fonction mathématique atteint sa valeur maximale ou minimale. En statistique, ce problème est souvent abordé en ces termes : trouver les valeurs des paramètres pour lesquels une fonction atteint son maximum ou son minimum.

Parfois, il est possible de trouver une solution algébrique à ce problème à l’aide du calcul différentiel et intégral. Par contre, il arrive qu’il soit trop difficile, voire impossible, de dériver la fonction en question. L’optimisation numérique est une bonne solution dans un tel cas.

Fonctions R utiles en optimisation numérique :

• pour optimiser une fonction à une variable : optimize,
• pour optimiser une fonction avec un nombre de variables quelconque : nlm, optim,
• optimisation sous contrainte : constrOptim.

Exemple d’optimisation d’une fonction à une variable en R : trouvons en quel point la fonction de densité de la loi normale atteint son maximum. La théorie nous dit que ce maximum est atteint en la valeur de l’espérance de la loi. Voyons si l’optimisation numérique saura retourner le bon résultat.

curve(
  expr = dnorm(x, mean = 3), from = -3, to = 10, 
  main = "Maximum de la densité normale", ylab = "densité N(3,1)"
)
abline(v = 3, lty = 2, col = "blue")

optimize(f = dnorm, interval = c(-3, 10), mean = 3, maximum = TRUE)

## $maximum
## [1] 3
## 
## $objective
## [1] 0.3989423

Oui, pour une loi normale d’espérance 3 et de variance 1, nous arrivons bien numériquement au résultat que le maximum de la densité est atteint en la valeur 3. La fonction optimize nous dit aussi que ce maximum vaut :

dnorm(3, mean = 3)

## [1] 0.3989423

Les fonctions nlm, optim et constrOptim utilisent des algorithmes itératifs. Elles ont besoin de valeurs initiales pour les paramètres (argument par à fournir obligatoirement). À chaque itération de l’algorithme, elles modifient ces valeurs en tentant de se diriger vers l’optimum de la fonction. Elles peuvent :

• ne pas converger,
• converger au mauvais endroit (optimum local plutôt que global).

Il faut être prudent lors de leur utilisation. Par exemple, optim est sensible au choix de plusieurs arguments, notamment :

• l’algorithme employé,
• les valeurs initiales données aux paramètres.

Ces fonctions sont tout de même très utiles pour effectuer une optimisation lorsque celle-ci est difficile ou impossible à réaliser algébriquement.

Voici un exemple d’optimisation d’une fonction à plusieurs variables. La fonction lm minimise le critère des moindres carrés, en implémentant des formules algébriques. Les estimations des paramètres du modèle linéaire que lm retourne sont les points en lesquelles la fonction des moindres carrés est minimisée. Tentons de minimiser cette fonction de façon numérique. Pour ce faire, nous avons d’abord besoin d’une fonction qui calcule le critère des moindres carrés et qui prend comme premier argument les paramètres du modèle. Nous n’avons pas encore vu dans le cours comment créer des fonctions, mais je me permets tout de même ici d’en créer une, pour illustrer l’optimisation numérique. La syntaxe pour créer des fonctions R sera vue au prochain cours.

Le critère des moindres carrés est calculé en sommant les différences au carré entre les valeurs observées d’une variable et les valeurs prédites par le modèle. Pour un modèle de régression linéaire, la fonction suivante calcule de façon matricielle la valeur du critère.

moindresCarres <- function(beta, y, X) {
  as.vector(crossprod(y - X %*% matrix(beta, ncol = 1)))
}

Le vecteur y doit contenir les valeurs observées de la variable réponse et la matrice X est la matrice de design du modèle. Cette dernière contient les observations des variables explicatives pour les termes présents dans le modèle. Le vecteur y et la matrice X sont des composantes du modèle supposées connues ici. C’est le vecteur de paramètre beta que nous cherchons à estimer. Nous allons utiliser les données du jeu de données cars dans cet exemple.

Voyons d’abord le résultat obtenu avec la fonction lm pour un modèle quadratique.

reg <- lm(dist ~ speed + I(speed^2), data = cars)
coefficients(reg)

## (Intercept)       speed  I(speed^2) 
##   2.4701378   0.9132876   0.0999593

Pour retrouver ce résultat par optimisation numérique, il faut d’abord construire le vecteur y et la matrice X comme suit.

y <- cars$dist
X <- cbind(intercept = 1, cars$speed, cars$speed^2)

La fonction lm arrive à la valeur minimale des moindres carrés suivante

moindresCarres(beta = coefficients(reg), y = y, X = X)

## [1] 10824.72

pour les valeurs de paramètres β = (2.4701378, 0.9132876, 0.0999593). Quel résultat est obtenu avec optim?

op1 <- optim(par = c(3,3,3), fn = moindresCarres, y = y, X = X)
op1

## $par
## [1] 7.2212674 0.2859028 0.1191485
## 
## $value
## [1] 10848.71
## 
## $counts
## function gradient 
##      144       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

L’algorithme a convergé (car il retourne une valeur de 0 pour l’élément convergence dans la sortie), mais il n’arrive pas au bon résultat.

Solution potentielle : changer d’algorithme d’optimisation.

op2 <- optim(par = c(3,3,3), fn = moindresCarres, y = y, X = X, method = "BFGS")
op2

## $par
## [1] 2.47011519 0.91329056 0.09995889
## 
## $value
## [1] 10824.72
## 
## $counts
## function gradient 
##       43        6 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Autre solution potentielle : changer les bornes initiales.

op3 <- optim(par = c(2.5,1,0.1), fn = moindresCarres, y = y, X = X)
op3

## $par
## [1] 2.46514183 0.91414522 0.09993205
## 
## $value
## [1] 10824.72
## 
## $counts
## function gradient 
##      150       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Ici, même en partant de valeurs initiales très proches des paramètres optimaux, l’algorithme d’optimisation utilisé par défaut avec optim n’arrive pas tout à fait à trouver l’optimum global de la fonction. Seule la solution de changer l’algorithme d’optimisation nous permet d’arriver approximativement au même résultat que celui trouvé algébriquement par lm.

coefficients(reg)

## (Intercept)       speed  I(speed^2) 
##   2.4701378   0.9132876   0.0999593

op2$par

## [1] 2.47011519 0.91329056 0.09995889

Cet exemple illustre comment la fonction optim s’emploie. Il faut d’abord lui donner en entrée des valeurs initiales pour les paramètres de la fonction à optimiser (argument par). Ensuite, il faut lui fournir la fonction R qui implémente la fonction mathématique à optimiser (argument fn). Cette fonction doit retourner une seule valeur, numérique. De plus, son premier argument doit obligatoirement être le vecteur des paramètres que nous cherchons à estimer par l’optimisation effectuée. Après les arguments par et fn, il faut fournir, au besoin, des arguments à passer à la fonction donnée en entrée via l’argument fn (les arguments y et X dans l’exemple). Finalement, nous pouvons configurer le fonctionnement de la fonction optim en modifiant les valeurs des arguments method, lower, upper, control, ou hessian.